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Hard spheres are ubiquitous in condensed matter: they have been used as models for liquids, crys- 
tals, colloidal systems, granular systems, and powders. Packings of hard spheres are of even wider 
interest, as they are related to important problems in information theory, such as digitalization 
of signals, error correcting codes, and optimization problems. In three dimensions the densest 
packing of identical hard spheres has been proven to be the FCC lattice, and it is conjectured that 
the closest packing is ordered (a regular lattice, e.g, a crystal) in low enough dimension. Still, 
amorphous packings have attracted a lot of interest, because for polydisperse colloids and granular 
materials the crystalline state is not obtained in experiments for kinetic reasons. We review here a 
theory of amorphous packings, and more generally glassy states, of hard spheres that is based on 
the replica method: this theory gives predictions on the structure and thermodynamics of these 
states. In dimensions between two and six these predictions can be successfully compared with 
numerical simulations. We will also discuss the limit of large dimension where an exact solution 
is possible. 

Some of the results we present here have been already published, but others are original: in 
particular we improved the discussion of the large dimension limit and we obtained new results on 
the correlation function and the contact force distribution in three dimensions. We also try here 
to clarify the main assumptions that are beyond our theory and in particular the relation between 
our static computation and the dynamical procedures used to construct amorphous packings. 
There remain many weak points in our theory that should be better investigated. We hope that 
this paper can be useful to present the state of the art of the method and to stimulate new research 
in this field. 
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I. INTRODUCTION 

The study of amorphous states of hard spheres is rel- 
evant for a large class of physical systems, including liq- 
uids, glasses, colloidal dispersions, granular matter, pow- 
ders, porous media. Therefore, after pioneering works 
done in the sixties (Bernal and Mason, 1960; Bernal 
et at, 1962; Finney, 1970; Mason and Clark, 1965, 1966; 
Scott, 1962; Scott and Kilgour, 1969), a huge amount 
of precise numerical and experimental data is now avail- 
able, see e.g. (Angelani and Foffi, 2007; Bennett, 1972; 
Berryman, 1983; Clarke and Jonsson, 1993; Donev et at, 
2005a; Jerkins et at, 2008; Krauth, 2006; Lubachevsky 
and Stillinger, 1990; Majmudar et at, 2007; Mathcson, 
1974; O'Hcrn et ai, 2002, 2003; Powell, 1979; Pusey 
and Van Megen, 1986; Rintoul and Torquato, 1996; Sil- 
bert et al., 2006; Skoge et al, 2006; Somfai et at, 2007; 
Speedy, 1998; Torquato, 1995, 2002). Moreover, the 



sphere packing problem is related to many mathematical 
problems and arises in the context of signal digitalization 
and of error correcting codes, and it has been investigated 
in detail by the information theory community (Conway 
and Sloane, 1993; Rogers, 1964). Nevertheless, a satisfac- 
tory characterization of the amorphous states of a system 
of identical hard spheres is not yet available and the def- 
inition of amorphous close packed states is still matter 
of debate (Kamien and Liu, 2007; O'Hcrn et at, 2002; 
Torquato et al., 2000). From the rigorous point of view, 
for space dimension d > 3 only some not very restric- 
tive bounds have been obtained, and in particular it is 
still unclear whether the densest packings for d — > oo arc 
amorphous or crystalline (see (Sloane, 2007) for a list of 
all known densest packings up to d = 128). 

Dense amorphous packings of hard spheres are usually 
produced according to some specific dynamical protocol. 
Typically one starts from an initial random configura- 
tion of the spheres, obtained e.g. by throwing them into 
a container, and then shake, tap, or agitate in some way 
the spheres until a jammed structure is found (Abate 
and Durian, 2006; Bennett, 1972; Daniels and Behringer, 
2006; Dauchot et at, 2005; Jerkins et ai, 2008; Majmu- 
dar et ai, 2007; Pica Ciamarra et at, 2007; Pusey and 
Van Mcgcn, 1986; Schroter et at, 2005; Scott and Kil- 
gour, 1969; Torquato, 2002). In numerical simulations, 
amorphous packings are produced by inflating the par- 
ticles while avoiding superposition either by molecular 
dynamics (Doncv et at, 2005a; Lubachevsky and Still- 
inger, 1990; Skoge et al., 2006) or by using soft particles 
and minimizing the energy (Clarke and Jonsson, 1993; 
O'Hcrn et ai, 2002, 2003; Silbcrt et at, 2006; Somfai 
et at, 2007). As a matter of fact, most of these proce- 
dures, if crystallization is avoided, lead to a final packing 
fraction close to 0.64 in d = 3 and to 0.84 in d = 2. These 
values of density, that are approximately 10% smaller 
than the values of the ordered close packing, have been 
called "random close packing density" . 

Unfortunately, the algorithms (or procedures) that are 
used to create such packings are complicated dynamical 
non- equilibrium procedures. Obtaining analytical results 
for the properties of the final states requires an analytical 
solution of such complicated dynamical processes, that 
is very difficult even in the simplest theoretical mod- 
els (Krzakala and Kurchan, 2007; Torquato and Still- 
inger, 2006a). The aim of this paper is then to identify 
a class of amorphous packings that might be described 
using equilibrium statistical mechanics, that is, in a static 
framework. These packings will be defined as the infinite 
pressure limit of glassy states of hard spheres: such glassy 
states, if dense enough, are well defined metastable states 
with very long life times, and should be then correctly 
described by equilibrium statistical mechanics. The idea 
of studying amorphous packings as the infinite pressure 
limit of a metastable state has been already discussed 
in the literature (Aste and Coniglio, 2004; Biroli and 
Mezard, 2001; Kamien and Liu, 2007; Krzakala and Kur- 
chan, 2007; Parisi and Zamponi, 2005; Pica Ciamarra 
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et al, 2003; Rivoire et al, 2004; Tarzia et al, 2004; Zam- 
poni, 2007) and is appealing because it converts a diffi- 
cult dynamical problem into a much simpler equilibrium 
problem. 

Our approach to study glassy states will be based on 
the so-called Random First Order Transition (RFOT) 
theory of glasses, whose theoretical foundations were 
posed in a series of papers by Kirkpatrick, Thirumalai 
and Wolynes (Kirkpatrick and Thirumalai, 1987; Kirk- 
patrick and Wolynes, 1987a), see (Cavagna, 2009) for 
a detailed and recent review. In this theory the glass 
transition of particle systems is assumed to be in the 
universality class of the 1-step Replica Symmetry Break- 
ing (1RSB) transition that happens in some mean-field 
exactly solvable spin glass models (Gross and Mezard, 
1984; Mezard et al, 1987). Under this assumption, the 
glassy states of realistic finite-dimensional systems can 
be studied analytically, within some approximation, us- 
ing equilibrium statistical mechanics by means of density 
functional theory (Chaudhuri et at, 2005; Dasgupta and 
Vails, 1999; Kim and Munakata, 2003; Kirkpatrick and 
Wolynes, 1987a; Singh et al, 1985; Stoessel and Wolynes, 
1984; Yoshidome et al., 2007) and of the replica trick 
(Mezard and Parisi, 1996; Monasson, 1995). In particu- 
lar, the replica method seems to give good quantitative 
estimates of the glass transition temperature (or density) 
and of the equation of state of the glass for Lennard- 
Jones systems (Coluzzi et at, 1999; Mezard and Parisi, 
1999a, b, 2000) and hard spheres (Cardenas et at, 1998, 
1999; Parisi and Zamponi, 2005, 2006a). 

Beside the large amount of numerical and experimen- 
tal data available in the literature, there is an important 
advantage in working with hard spheres with respect to 
Lennard-Jones-like potentials. "Ground states" of hard 
spheres are obtained in the infinite pressure limit, and 
correspond then to sphere packings that have interesting 
geometrical and topological properties, that can be in- 
vestigated by looking at the network of contacts and at 
contact forces. This is very interesting because one can 
identify a set of geometrical observables that can be com- 
puted within the theory and directly compared with sim- 
ulation and experiments. In this paper we will show that 
this allows a very precise test of the theory. Moreover, a 
more direct geometrical approach is possible; it has been 
largely exploited, sec e.g. (Anikccnko and Mcdvcdcv, 
2007; Anikeenko et at, 2008; Aste, 2005; Aste et al, 2005; 
Donev et al., 2007, 2005a; Kansal et al., 2002; Lochmann 
et al, 2006; van Meel et al, 2009; Richard et al, 1999; 
Skoge et al, 2006), and led to some important successes 
in the characterization of sphere packings. This is not 
possible for Lennard-Jones-like particles since the (zero 
temperature) ground states of the system do not have 
special geometrical properties, and indeed their struc- 
ture is quite similar to typical liquid configurations. One 
might consider instead a soft potential that vanishes out- 
side a finite radius; then the zero-energy ground states 
correspond to hard sphere configurations. In this way 
one obtains a system that displays the same geometri- 



cal properties of hard spheres at zero temperature and 
energy, but at the same time becomes soft at finite tem- 
perature. This has been largely exploited (Berthier and 
Witten, 2009a,b; O'Hern et al, 2002, 2003; Schrcck and 
O'Hern, 2008) to obtain important informations on the 
properties of amorphous hard sphere packings. Therefore 
in this paper we will focus on the hard sphere case and we 
will discuss in details the limits of the theory and how it 
compares with numerical and experimental results. We 
will show that despite the strong idealizations involved 
in the theory, the agreement with numerical data is sur- 
prisingly good. Moreover, we will be able to study in full 
detail the limit of large space dimension, and obtain the 
asymptotic value of the density of amorphous packings. 

Remarkably, a class of mean-field hard sphere mod- 
els have been recently formulated, for which the RFOT 
scenario is exact (Biroli and Mezard, 2001; Krzakala 
et al, 2008; Mari et al, 2009; Pica Ciamarra et al, 2003; 
Rivoire et al, 2004; Sellitto et al, 2005; Tarzia et al, 
2004). These models allowed to test the methods used 
here confirming that they are reliable, at least at the 
mean- field level. In particular, (Mari et al, 2009) for- 
mulated a model in this class whose phase diagram is 
exactly the same as the one we will discuss below for 
finite dimensional hard spheres. 

It is worth to stress that in experiments on granular 
systems and powders the role of friction is very impor- 
tant (Abate and Durian, 2006; Daniels and Behringer, 
2006; Dauchot et al, 2005; Lechenault et al, 2008; 
Pica Ciamarra et al, 2007; Schrdter et al, 2005; 
Shundyak et al, 2007; Somfai et al, 2007), for instance 
in determining the existence of loose packings (Jerkins 
et al, 2008; Onoda and Liniger, 1990; Song et al, 2008). 
Friction complicates a lot the theoretical analysis of the 
packing problem, since the system is intrinsically out of 
equilibrium and standard equilibrium statistical mechan- 
ics is in principle useless. Nevertheless, since the pioneer- 
ing work of Edwards (Edwards and Oakeshott, 1989; Ed- 
wards, 1998), statistical mechanics ideas have been used 
to describe frictional systems, leading to remarkable re- 
sults (Goldbart et al, 2005). The comparison with ex- 
perimental results is made difficult by the fact that in 
most experiments samples are polydispcrse, often with a 
quite large range of particle sizes as in the case of many 
granulars. 

For reasons of space, this paper will be focused on our 
approach that we wish to discuss in full detail; therefore 
in the following we will consider mainly the statistical 
properties of a system of frictionless spheres since our 
method is based on equilibrium statistical mechanics. We 
will not discuss in detail neither the geometrical proper- 
ties (unless needed to compare numerical data with our 
results) of amorphous packings, nor how their properties 
are influenced by the presence of friction. While in prin- 
ciple polydispersity can be included in our theory, here 
we limit ourselves to discuss the simple case of monodis- 
perse systems, and only at the end we consider the case 
of binary mixtures. The literature on hard spheres in 
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general (and on the role of friction and geometry in par- 
ticular) is immense and covering it here would require a 
considerable effort: the reader is therefore referred to the 
original literature and to the existing excellent books and 
reviews, e.g. (Alexander, 1998; Aste, 2005; Conway and 
Sloane, 1993; Goldbart et at, 2005; Krauth, 2006; Liu 
and Nagel, 2001; Rogers, 1964; Torquato, 2002). Simi- 
larly, many excellent general books and reviews on the 
physics of glasses and the glass transition exist, see e.g. 
(Binder and Kob, 2005; Cavagna, 2009; Debenedetti and 
Stillinger, 2001; Donth, 2001; Ediger et at, 1996; Leuzzi 
and Nieuwcnhuizcn, 2007) just to quote a few of them. 



A. Organization of the paper (and how to read it) 

This paper is organized as follows. Section II is de- 
voted to a general discussion of the ideas that lead to 
the connection between packings that are produced by 
dynamic protocols and infinite pressure glassy states. In 
section III the replica method is introduced in a general 
context and it is explained how it can be used to com- 
pute the phase diagram of glassy states. In section IV a 
first implementation of the method, that works far from 
jamming, is discussed. Sections V, VI, VII, VIII contain 
the core of the paper: the method is implemented in a 
way that works up to jamming and most of the results 
are presented. In particular section VI contains the dis- 
cussion of the d — > oo limit. Many technical parts of the 
paper are in the Appendix. 

We tried where possible to make the different parts 
of the paper independent. There arc different ways, of 
increasing difficulty, to read the paper: 

• The discussion of section II is self-contained. We 
suggest to the reader not interested in technical de- 
tails to read only section II, then skip the compu- 
tations and look directly to the figures and tables 
that contain most of our results, and finally jump 
to the conclusions. Each figure where results are 
presented contains a reference to the section where 
the corresponding calculations are discussed. 

• The general discussion of the replica method in sec- 
tion III is as self-contained as possible. This might 
be read together with Appendix A if one wants to 
obtain more insight into the method without going 
into the technical details of the computations. 

• The technical sections IV, V, VI, VII, VIII are all 
independent one from the other; to understand one 
of them in detail one needs only to read section III 
and the Appendices. 

• Finally, in IX we present some results on the ex- 
tension of the theory to binary mixtures, that are 
based on the method of section VII. 

It is important to stress that we decided to present 
different implementations of the method separately, to- 
gether with the corresponding results. Clearly, another 



possible choice could be to collect all the results together 
in a separate section. Our choice has been made to stress 
that each approximation scheme gives slightly different 
results; although the qualitative picture stays the same, 
it is difficult to compare quantitatively different approxi- 
mations. Moreover, each approximation has some advan- 
tages and disadvantages that we tried to discuss in detail; 
in particular, some observables can be computed within 
some approximation and not within others. We leave to 
the reader the task to compare the different methods, and 
choosing the best one according to his personal taste. 



B. Notations 

It is useful to give here some definitions that will be 
widely used in the following. We will consider a system 
of hard hypcrspheres in d dimension with diameter D. 
We will denote by 



V d (D) 



T d/2 



r(i + d/2) 



D" 



^(D)=dD^V d (D) = ^D^ 



(1) 



respectively the volume and surface of a sphere of radius 
D. It is also convenient to define 



V d = V d (l) 



T d/2 



r(i + d/2) 

27T d / 2 

r(d/2) ' 



(2) 



Note that fl d is the d-dimensional solid angle. 

Often, if there is no ambiguity, we will use the symbols 
x, y, ■ ■ ■ to denote vectors in R d . If there is an ambigu- 
ity, we will use the letter r and in this case we will use 
r for a <i-dimensional vector and r = \r\ for its modulus. 
Correspondingly / dfF(r) will indicate an integral over 

r € JR d while drF(r) will indicate an integral over the 
real r. We will use the notation S(\r\ — D) for the dis- 
tribution defined by / drF(\r\)S(\r\ - D) = T, d (D)F(D), 
while 6(r-D) is defined by f °° drF(r)5(r-D) = F(D). 
They are related by 5(\r\ - D) = n d S(r - D). 

For a generic potential <j)(r) we define b(r) = e~^ r ) 
and the Mayer function f(r) = b(r) — 1. For the specific 
case of Hard Spheres we have b(r) = x(r) = 8(\r\ — D). 

Defining p = N/V the number density of the spheres, 
we introduce as usual the packing fraction ip = pV d (D/2), 
i.e. the fraction of volume covered by the spheres. In the 
following, when talking about "density" , we will usually 
refer to the packing fraction. 



II. A CLASS OF AMORPHOUS PACKINGS: INFINITE 
PRESSURE GLASSY STATES 

In this section we will define a class of "amorphous 
packings" of hard spheres that can be studied within the 
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FIG. 1 (From (Skoge et ai, 2006)) Evolution of the pressure during compression at rate 7 in d = 3 (left) and d — 4 (right). 
The density ip is increased at rate 7 and the reduced pressure p(p) = fiP/ p is measured during the process. See (Skoge et al., 
2006) for details. The quantity (p~j{(p) = ^g^~3 is plotted as a function of tp. If the system jams at density ipj, p — ¥ 00 and 
(fij — > ipj. Thus the final jamming density is the point where Pj{p) intersects the dot-dashed line tpj — p. (Left) The dotted 
line is the liquid (Percus-Yevick) equation of state. The curves at high 7 follow the liquid branch at low density; when they 
leave it, the pressure increases faster and diverges at <pj. The curves for lower 7 show first a drop in the pressure, which signals 
crystallization. (Right) All the curves follow the liquid equation of state (obtained from Eq.(9) of (Bishop and Whitlock, 2005)) 
and leave it at a density that depends on 7. In this case no crystallization is observed. For 7 = 10 -5 the dot-dashed line is a 
fit to the high-density part of the pressure (glass branch). The arrow marks the region where the pressure crosses over from 
the liquid to the glass branch. 



framework of equilibrium statistical mechanics. We start 
by discussing some algorithms that are commonly used to 
construct amorphous packings. Then we argue that the 
final states reached by these algorithms are well defined 
metastable states whose properties can be investigated 
by a "static" computation (i.e. without any knowledge of 
the dynamical process that generated the packings) . This 
point is very delicate and is receiving a lot of attention 
in the context of optimization problems (Krzakala and 
Kurchan, 2007), where it has not yet been solved. The 
discussion that follows is tentative and the problem surely 
deserves further investigation. 



A. Algorithms to construct amorphous packings 

The usual way to construct amorphous packings in ex- 
periments or numerical simulations is to compress the 
system according to some protocol. In early experiments 
particles were thrown randomly in a box which was then 
shaken (Scott and Kilgour, 1969), or were deposited ran- 
domly around a small seed cluster (Bennett, 1972). In 
numerical simulations a common protocol (Lubachevsky 
and Stillinger, 1990) is to slowly increase, at a given rate 
7, the particle diameter during a molecular dynamics 
run; it has been recently used extensively in three (Doncv 
et ai, 2005a) and higher (Skoge et ai, 2006) dimensions 
to produce amorphous packings of N ~ 10 4 spheres. An- 
other possibility is to increase the diameter of the spheres 
until two of them overlap, then eliminate the overlap by 



following a gradient descent using some potential vanish- 
ing outside the radius of the particle (Gao et at, 2006; 
Xu et al, 2005); or, alternatively, to start from a random 
configuration and minimize the energy at fixed density, 
repeating the procedure while increasing the density un- 
til it becomes impossible to find a zero-energy final con- 
figuration (O'Hern et ai, 2002, 2003). These two pro- 
cedures give very similar results (Schreck and O'Hern, 
2008). Other similar algorithms have been proposed and 
analyzed in (Clarke and Jonsson, 1993; Jodrey and Tory, 
1985; Lochmann et al, 2006). 

Based on standard concentration arguments, it is be- 
lieved that, in the limit N — > 00, the density of the final 
state is independent on the randomness built in the algo- 
rithm 1 (e.g. the initial configuration). This has been nu- 
merically verified for the soft-potential algorithms (Gao 
et ai, 2006; O'Hern et al, 2002, 2003; Xu et ai, 2005); 
the final density is very close to 0.64 and has been called 
J-point. It is a fact that for all algorithms that have been 
devised to construct amorphous packings of monodis- 
perse frictionless spheres, for N —¥ 00 the final density 
of the system converges to a value of <p close to 0.64 in 
three dimensions. It has been proposed to call the value 
(p = 0.64 random close packing density and different more 
precise definitions of this concept have been proposed in 



Similar results have been shown for some classes of algorithms 
to solve optimization problems, see e.g. (Barthel et al, 2003; 
Semerjian and Monasson, 2003). 
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the literature. However the precise numerical value of 
the latter quantity is found to depend on the details of 
the experimental protocol, and this led some authors to 
criticize the notion of random close packing (Torquato 
et al, 2000). 

To illustrate this difficulty, let us concentrate on data 
obtained using the Lubachevsky-Stillinger procedure at 
different rates 7 (Skoge et at, 2006) and reproduced in 
figure f . In this algorithm, during the compression, the 
(hard) spheres evolve according to molecular dynamics at 
some temperature; the value of the temperature is irrele- 
vant and fixes the unit of time (Skoge et al, 2006). One 
can measure the (kinetic) reduced pressure p(<p) = (3P/p, 
where /3 = 1/T, during the evolution. This quantity is 
reported in figure I for different values of 7. One observes 
that, on compressing the low density liquid at a constant 
rate 7, the pressure of the system follows the equilibrium 
pressure of the liquid up to some density (p g {~f) (often 
called glass transition density) around which the pressure 
starts to increase faster than in equilibrium, and diverges 
on approaching a value of density ^(7), which is called 
jamming density. At this point the algorithm stops be- 
cause the system cannot be compressed anymore: most 
of the spheres are in contact with their neighbors 2 . Val- 
ues of ^(7) have been accurately measured in (Skoge 
et at, 2006) as a function of 7. On the contrary, <p g (~i) 
is not precisely determined as long as 7 > 0: the glass 
transition is smeared and happens in a crossover region 
[ip~ (7), ^"(7)] (marked by an arrow in the right panel 
of figure 1). However, the amplitude of the crossover 
interval seems to decrease 3 for 7 — > 0, see figure 1. 

To characterize the typical configurations reached by 
the algorithm at (p = ^(7), one could try to solve 
the dynamics of the algorithm and compute for instance 
the correlation function g(r) as a function of "time". 
This has been done for much simpler algorithms, such 
as the Ghost Random Sequential Addition ( GRSA ) algo- 
rithm (Torquato and Stillinger, 2006a), which however 
are unable to reach interesting densities 4 . More effi- 
cient algorithms such as the Random Sequential Addi- 
tion (RSA), where one attempts to add a sphere ran- 
domly and accepts the move only if there are no over- 
laps, already cannot be analyzed analytically and one 
has to resort to numerical investigation (Talbot et at, 
2000; Torquato et at, 2006). 

For this reason, in order to compute analytically the 
properties of jammed configurations, one would like to 
make use of a "static" computation. This can be justi- 
fied as follows: if we plot the jamming density <Pj(j) as a 



2 At low compression rates, crystallization is observed in d = 3, but 
it seems to be strongly suppressed by kinetic effects in dimension 
d > 3 so we will neglect it for the moment. 

3 See (Mollcr et al., 2006) for a recent theoretical discussion of 
these effects. 

4 For instance in d = 3 the limiting density for the algorithm is 
ip = 0.125. 



function of 7, we obtain the plot reported in right panel 
of figure 2. For 7—5-0, the algorithm is equivalent to 
an equilibrium compression and the pressure should fol- 
low the equilibrium equation of state; then the final state 
will be the most dense state, which is a crystal at least if 
the dimension is not very large. However, for small but 
non-zero 7, crystallization is not observed and the data 
of (Skoge et al., 2006) suggest the existence of a plateau 
at some value of density, that we will call <pgcp for rea- 
sons that will be clear in the following. This is a hint 
of the existence of a long-lived metastable state; if this 
is the case, one can compute its properties by mean of a 
thermodynamic computation by restricting the partition 
function to the region of amorphous configurations. In 
section II. B we discuss this issue in more detail. 

It seems that many different algorithms that produce 
packings starting from a random configuration and with- 
out allowing the particles to relax much (Gao et al, 
2006; O'Hern et al., 2002, 2003; Xu et al, 2005) lead 
to very close values of density (the J-point). There- 
fore for large 7 (but still small enough to allow for 
some minimal relaxation of the particles) the final den- 
sity of the Lubachevsky-Stillinger algorithm should be 
close to the one reached by the J-point procedure, lead- 
ing to the schematic plot reported in left panel of fig- 
ure 2. This is what indeed happens in a simple mean- 
field model (Cugliandolo and Kurchan, 1993), and seems 
confirmed by the data for spheres in d — 4 reported in 
the right panel of figure 2. However, the validity of this 
statement is very debated for more general mean-field 
models (Krzakala and Kurchan, 2007). We will come 
back to this point in the following. 

Before concluding this section, let us stress that here 
we focused only on the behavior of algorithms that are 
currently used to construct amorphous packings; these 
algorithms typically perform local moves in phase space. 
Many other algorithms to simulate hard spheres have 
been invented. Smarter algorithms can be designed, that 
are able to sample the equilibrium measure at higher den- 
sity (Dress and Krauth, 1995; Grigera and Parisi, 2001; 
Krauth, 2006). 



B. The "equilibrium" phase diagram 

The equilibrium phase diagram of hard spheres in M 3 
is sketched in Fig. 3, where the pressure is reported as a 
function of the packing fraction. At low density the sys- 
tem is in a liquid phase (as defined e.g. by the low density 
virial expansion); the maximum possible density is real- 
ized by the FCC lattice, as conjectured by Kepler and 
rigorously proven by Hales (Hales, 2005). A first order 
phase transition between the liquid phase and the FCC 
crystal phase is found by numerical simulations (Alder 
and Wainwright, 1957; Wood and Jacobson, 1957) and 
in experiments on colloidal systems (Phan et al, 1996; 
Pusey and Van Megen, 1986). This is the equation of 
state that the system will follow if compressed at really 
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FIG. 2 (Left) A schematic plot of the jamming density v?j(7) as a function of compression rate 7 for the Lubachevsky-Stillinger 
algorithm. At quite high rates, the algorithm should be similar to the J-point procedure and converge to packings of density 
ipjif large) ~ <pj. On decreasing 7 the jamming density decreases (compare with figure 1) and reaches a plateau at ipj — ipgcp- 
For smaller 7 the system is able to crystallize and <Pj(0) = tp C rystal- In d = 3, crystallization is fast and the plateau is not 
observed, while in d > 3 crystallization is so slow that it is not observed at all in numerical simulations (Skoge et al., 2006). 
(Right) The same plot for hyperspheres in d = 4. The full line is the close packing density achieved by the D4 crystal, 
V?D4 = 7r 2 /16 (Conway and Sloane, 1993). Circles are values of </'j(7) (the point where the curves intersect the dashed line 
in right panel of figure 1). The dot-dashed line is the value of the J-point density (Schreck and O'Hern, 2008). In addition, 
we report the estimates of the glass transition density <p g ("f) (squares); error bars mark the amplitude of the crossover region 
and correspond to the arrow in right panel of figure 1. Dashed lines are fits to (Pj("y) = 0.473 + 0.023/(log 10 (7) + 0.85) and 
<p B (i) = 0.409 + 0.02/(log 10 ( 7 ) + 2.13). 



infinitesimal rate. We wish instead to focus on a small 
but finite rate in such a way to follow the liquid branch 
of the equation of state inside its metastability region. 



1. What is the fate of the liquid above freezing? 

The first idea to define amorphous states of hard 
spheres at high density is to assume that the liquid phase 
can be continued above the freezing density ipf and to 
look at its properties at large density (Aste and Coniglio, 
2004; Kamien and Liu, 2007). For instance one can 
choose a functional form that represents well the equa- 
tion of state of the liquid below ipf (e.g. the Carnahan- 
Starling or Percus-Yevick equation of state (Hansen and 
McDonald, 1986)) and assume that it describes the liq- 
uid phase also above <pf. On increasing the density the 
pressure of the liquid increases as the average distance 
between particles decreases: one may expect that it di- 
verges at a point where the particles get in contact with 
their neighbors, and the system cannot be further com- 
pressed. Then one might identify this point with the 
random close packing density, see left panel of Fig. 3. 

The first objection that has been raised against this 
proposal is that an intrinsic stability limit of the liquid (a 
spinodal point) might exist at a density above ipf due to 
thermodynamic or kinetic reasons. It is probable that a 
thermodynamic spinodal does not exist, because any rea- 
sonable continuation of the liquid equation of state does 
not predict such an instability (manifested e.g. by an in- 



finite compressibility). A kinetic spinodal, related to the 
existence of the crystal (Cavagna et al, 2005), could in- 
stead exist, at least in monodisperse systems. This would 
imply the impossibility of reaching amorphous jammed 
states if the compression rate is not very high. We will 
assume in the following that the metastable liquid can 
be compressed as slow as wished avoiding crystallization. 
This is not a very good assumption for monodisperse 
three-dimensional spheres but seems to be more close to 
reality for larger dimension (Skoge et al., 2006), see fig- 
ure 1, or for suitably chosen binary mixtures. This point 
requires a better investigation, e.g. following the analy- 
sis of (Cavagna et al, 2005), and we will not discuss it 
further here. 

A second objection is that in presence of a first or- 
der phase transition the continuation of one phase into 
its metastable region is not well defined due to the ap- 
pearance of essential singularities at (pf. Many possible 
continuations of the low density equation of state above 
iff are possible: the properties of the "liquid" above tpf 
depend on the history of the sample, much as it happens 
for the hysteresis of a ferromagnetic system. Neverthe- 
less, the ambiguity is expected to be exponentially small 
in the distance from tpf, and as the distance between 
tpf and the maximum density lpfcc is not so large, one 
might expect to obtain a meaningful result in this way. 
And indeed different possible continuation of the liquid 
equation of state (e.g. Carnahan-Starling, Percus-Yevick, 
Hypernetted Chain) differ by less than 10% in the dense 
region around <p ~ 0.64 in d — 3. Moreover, the am- 
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FIG. 3 Schematic phase diagram of hard spheres in M . (Left) Continuation of the liquid equation of state in the metastable 
region. (Right) Expected behavior in presence of a thermodynamic glass transition. 



biguity becomes smaller and smaller on increasing the 
dimension. 

Still the extrapolation to infinite pressure can give very 
different results; in (Kamien and Liu, 2007) it is shown 
that one can reasonably fit the pressure in the metastable 
region by a free volume equation of state, and obtain a 
divergence of the pressure at ifincp ~ 0.64 in d = 3. 
On the other hand, all possible analytic continuations 
of the liquid equation of state, based on resummations 
of the virial series, predict a divergence of the pressure 
at unphysical large values of ip: e.g. the CS equation 
predicts a divergence in <p = 1 which is clearly wrong 
since it is larger than the FCC value and moreover im- 
plies that the available volume is completely covered by 
the spheres. Moreover, the g(r) computed using these 
resummations do not show the characteristic features ob- 
served in amorphous jammed packings. Therefore, from 
a theoretical point of view it is very difficult to justify the 
validity of the phenomenological free volume fit discussed 
in (Kamien and Liu, 2007) on the basis of an analytical 
continuation of the liquid equation of state. Clearly we 
cannot exclude that a more refined resummation of the 
virial series will give a continuation of the liquid branch 
with a divergence in 0.64 and a g(r) similar to the one of 
amorphous packings; but such theory has not been found 
yet. This motivates the study of different proposals, such 
as the existence of a phase transition in the liquid branch. 
This idea is also supported by accurate integral equations 
based on improvements of the PY equation that seem to 
indicate the existence of a phase transition at high den- 
sity (Robles and de Haro, 2003). 



2. The ideal glass transition 

In this section we discuss the existence of a thermody- 
namic glass transition in the metastable liquid branch of 
the phase diagram (Cardenas et al, 1998; Speedy, 1994, 
1998; Stoessel and Wolynes, 1984; Woodcock and An- 



gell, 1981). Such a transition can be expected for vari- 
ous reasons: it is predicted by mean-field models (Biroli 
and Mezard, 2001; Pica Ciamarra et al, 2003; Rivoire 
et al, 2004), and is suggested by the results of (Rob- 
les and de Haro, 2003). Moreover, from a dynami- 
cal point of view, an ergodicity breaking transition is 
predicted by Mode-Coupling theory (Bengtzelius et al, 
1984; Gdtze, 1999; van Megen and Underwood, 1993; Sel- 
litto et al, 2005), and a recently published set of accu- 
rate experimental and numerical data (Berthier and Wit- 
ten, 2009a, b; Brambilla et al, 2009) on three dimensional 
spheres strongly indicates a divergence of the equilibrium 
relaxation time in the metastable liquid phase at a den- 
sity ipo at which the pressure is still finite. This sup- 
ports the existence of a phase transition at ipo in three 
dimensions. In two dimensions the situation seems quite 
different (Donev et al, 2006; Santen and Krauth, 2000; 
Tarzia, 2007). 

Still, the existence of an ideal glass transition in a 
finite-dimensional system is a matter of intense debate, 
see e.g. (Bouchaud and Biroli, 2004; Brumer and Re- 
ichman, 2004; Cavagna, 2009; Donev et al, 2006; San- 
ten and Krauth, 2000; Tarzia, 2007; Xia and Wolynes, 
2001a, b). To avoid entering in this delicate discussion, 
for the moment, we will take a more "pragmatic" point 
of view. We will assume that a thermodynamic glass 
transition exists and investigate its consequences. It will 
become clear in a moment that for the comparison with 
numerical and experimental data the existence of a true 
thermodynamic transition is not a key issue, since in real 
life one is always stuck in metastable glassy states. 

So let us assume again that the liquid phase can be 
continued above tp / and neglect the (small) ambiguity in 
its definition due to its intrinsic metastability with re- 
spect to crystallization. We assume that at a density ipx 
a thermodynamic glass transition (sometimes called ideal 
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glass transition) happens 5 . The transition is signaled by 
a jump in the compressibility of the system. A simple 
qualitative argument to explain this is the following: in 
the dense liquid phase particles vibrate on a fast time 
scale in the cages made by their neighbors, while on a 
much larger time large scale cooperative relaxation pro- 
cesses happen (structural relaxation). If we change the 
density by Aip the pressure will instantaneously increase 
by a APq ; very rapidly the average size of the cages will 
decrease a little due to the increase in density and the 
pressure will relax to a value APy < AP - Then, on 
the time scale of structural relaxation, the structure will 
change to follow the change in density and the pressure 
will relax further to a value APqo < APf. We assume 
that at the glass transition the latter relaxation is frozen 
and the corresponding time scale becomes infinite: in 
other words, as usually done in the glass transition litera- 
ture (Ediger et ai, 1996), we identify ipx with the density 
ipo defined in (Berthier and Witten, 2009a, b; Brambilla 
et ai, 2009). Thus, in the glass phase the increase in 
pressure following a change in density will be larger than 
in the liquid phase, leading to a smaller compressibility 
K = p-^Ap/AP). 



The schematic phase diagram that we expect in pres- 
ence of a glass transition is in right panel of Fig. 3. The 
existence of a glass transition can cure the behavior of 
the pressure of the liquid, that seems to diverge at a den- 
sity bigger than the FCC density. The pressure of the 
glass diverges at a "Glass Close Packing" density <pgcp 
and we could be tempted to identify ipccp with the RCP 
density. In next section we discuss a further complica- 
tion: the existence of a very large number of glassy states, 
in addition to the ideal (thermodynamic) glass, with dif- 
ferent densities. Before concluding this section, we wish 
to stress again that, although recent very accurate nu- 
merical data (Berthier and Witten, 2009b) support the 
existence of a glass transition of the type discussed in this 
section for Hard Spheres, this remains one of the most 
debated problems in the community. 



The density ipx has been associated with the name of Kauzmann, 
who first observed that in some molecular glasses, extrapolating 
the liquid entropy at low temperature, one would obtain at some 
point an entropy smaller then that of the crystal (Kauzmann, 
1948). Based on the intuitive notion that a liquid should have 
more states than a crystal, hence more entropy, he concluded 
that some kind of transition should happen preventing the liquid 
entropy to cross the crystal one. However, in the case of hard 
spheres, it is well known that the crystal entropy becomes bigger 
than that of the liquid at the freezing/melting transition. Hence 
the entropy of the liquid is smaller than that of the crystal in the 
whole metastable liquid branch. At (fx the complexity (that will 
be introduced below) vanishes; but this quantity has nothing to 
do with the difference between liquid and crystal entropies in the 
case of hard spheres. 



3. Many glassy states: the "mean-field" phase diagram 

The glass transition, in the standard picture coming 
from the analysis of mean field models, is related to the 
appearance of many metastable glassy states 6 in addi- 
tion to the ideal glass one. These states appear in the 
liquid above some density ipd and can be defined for in- 
stance as minima of a suitable density functional (Chaud- 
huri et ai, 2005; Dasgupta and Vails, 1999; Kirkpatrick 
and Wolynes, 1987a; Mezard et ai, 1987; Thoulcss et at, 
1977) (see Appendix A for a more detailed discussion). In 
the interval of densities ipd < ip < tfK, particles vibrate 
around these locally stable structures, that are visited 
subsequently on the scale of the structural relaxation. 
If we "artificially" freeze the structural relaxation and 
compress the system, the pressure increases faster than 
if the structure is allowed to relax for the same reason as 
above: the system is forced to reduce the size of the cages 
to respond to a change in density. The pressure diverges 
at the point where the particles get in contact with their 
neighbors and the average size of the cages is zero. In this 
way, to each configuration of the liquid at a given density 
<p € [ipd,<PK] one can associate a jammed configuration 
at a density tpj(tp) that is obtained by compressing this 
configuration fast enough to avoid structural relaxation. 
A "glass state" can be roughly thought as a set of config- 
urations leading to the same jammed configuration after 
a fast compression 7 (Speedy, 1998). 

In the Lubachevsky-Stillinger protocol discussed be- 
fore, (see Fig. 1) one chooses a priori a compression rate 
7 and the system can equilibrate only up to a density 
fgil) where the relaxation time becomes of the order of 
the compression rate. At this density the structure can- 
not follow anymore the compression and basically the sys- 
tem responds to the compression by reducing the cages 
up to the jamming density Calling ip t h the value 



There arc here two "types" of metastability: the first is the 
metastability of the whole amorphous branch with respect to 
crystallization. We are now discussing a second type of metasta- 
bility, i.e. the existence of glassy states that have lower density 
(or smaller entropy) with respect to the glassy state with maxi- 
mal density (the ideal one) and consequently they are metastable 
with respect to it. In order to avoid confusion, we will indicate 
that some of the glassy states are metastable with respect to the 
ideal glass only when needed. Otherwise, we will use "glassy 
states" to indicate both metastable glasses and the ideal glass. 
In any case these states will have a small, but finite probability 
of decaying into the ideal state; this probability should go ex- 
ponentially to zero when their density approaches the density of 
the ideal glassy state. 

In systems with soft potentials a similar procedure has been pro- 
posed in (Stillinger and Weber, 1982, 1985) in order to associate 
to each configuration of energy e a mechanically stable config- 
uration or "inherent structure" of lower energy ejg: starting 
from the reference configuration one quenches the system at zero 
temperature and finds a minimum of the potential, which is the 
corresponding inherent structure. We are describing the same 
procedure, with energy replaced by (inverse) density and tem- 
perature by (inverse) pressure. See (Speedy, 1998) for details. 
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of the jamming density of the states that hrst appear at 
ipd, it follows that (in the mean held picture) we can pro- 
duce jammed configurations in a whole range of densities 
<Pth < f < fGCP (Krzakala and Kurchan, 2007; Zam- 
poni, 2007). The notation ^gcp (Glass Close Packing 
density) is appropriate since <pgcp is the highest possi- 
ble density of infinite pressure glassy states. 

Note that, as the structural relaxation time scale is ex- 
pected to diverge on approaching ipK (Berthier and Wit- 
ten, 2009a, b; Brambilla et al, 2009), at some point it 
will fall beyond any experimentally accessible time scale 
and necessarily the system will be frozen into a glassy 
state which is not the ideal glass. The ideal glass states 
are unobservable in practice: one will observe instead 
a nonequilibrium glass transition to a state that, again, 
will depend on the experimental protocol. This is very 
familiar in the structural glass literature, and in fact the 
diagram in hgure 4 (left panel) is analogous to the usual 
specific volume vs temperature plot (right panel in fig- 
ure 4) that characterizes structural glasses, if one identi- 
fies ip with the inverse of specific volume and 1/P with 
temperature. 

To summarize, the phase diagram inspired by mean- 
field models is reported in Fig. 4 and in table II. It is 
characterized by the presence of the following phases and 
transitions: 

• At equilibrium, a low density-liquid phase and a 
high-density crystal phase separated by a first- 
order transition, with corresponding freezing den- 
sity ip f and melting density ip m (black full line and 
black dots in Fig. 4). 

• A metastable dense liquid phase; above some den- 
sity (fid this phase is made by a collection of glassy 
states corresponding to locally stable configura- 
tions around which the system vibrates for a long 
time. Thus, ipd is defined, in a static framework, as 
the density at which glassy states first appear 8 and 
for this reason is often called clustering transition 
density 9 . 

• The liquid exists up to a density ifiK where an ideal 
glass transition happens. The ideal glass transition 
is signaled by a jump in the compressibility and by 
a divergence of the equilibrium relaxation time of 
the liquid 10 (Berthier and Witten, 2009a; Brambilla 
et al, 2009). 



At the mean-field level this point corresponds to the Mode- 
Coupling transition <Pmct 

Glassy states are called clusters in optimization problems 
Mode-Coupling Theory would predict a divergence of the relax- 
ation time at the smaller density tfiMCT ~ fd\ however, in a 
finite dimensional system there arc very strong arguments indi- 
cating that the liquid can be equilibrated up to ipx, thanks to 
activated processes that allow to jump over the barriers separat- 
ing metastable glassy states, see Appendix A. 



• For each density ip € [fd, fx], & different group 
of glassy states dominate the partition function. 
These can be followed by compressing very fast, 
and each group is characterized by a jamming den- 
sity We call fGCP = <Pj(<PK) the jamming 
density of the ideal glass state, and tpth — fjifd) 
the jamming density of the less dense glassy states. 

The number of glassy states corresponding to jamming 
densities tpj G [fth, feep] grows exponentially with the 
size of the system, N{(fj) — exp[NY,(ipj)]. The func- 
tion S(Vj)j usually called complexity or configurational 
entropy 11 , has the shape reported in the inset of Fig. 4: 
it jumps at a positive value at tpth and vanishes at <Pgcp 
(therefore, the number of ideal glasses is not exponen- 
tial in the system size). A more precise definition of this 
equilibrium complexity is presented in Appendix A. 



4. Remarks on the static phase diagram 

The striking similarity of the phase diagram in Fig. 4 
with the numerical results 12 of (Skogc et at, 2006), see 
Fig. 1, makes it a good starting point to understand 
the physics of jammed amorphous states. However one 
should keep in mind the following important remarks: 

1 . This phase diagram is an idealization that discards 
many difficulties in the definition of amorphous 
packings, mainly related to metastability, either of 
metastable glass states with respect to the ideal 
glass state, or of the ensemble of glassy states with 
respect to the crystal. The ambiguity in the def- 
inition of the liquid equation of state due to its 
metastability (Kamien and Liu, 2007) affects also 
the glass: the glass equation of state is theoreti- 
cally not well defined, with an ambiguity of the or- 
der of 10%, depending on the equation of state one 
chooses to describe the liquid. Note that these dif- 
ficulties might be not so important in some cases 
where the nucleation time of the crystal is very 



11 Although sometimes used as synonyms, in other cases these 
words denote different concepts: in fact, especially in the ex- 
perimental literature the word "configurational entropy" denotes 
the difference between the entropy of a system in its liquid and 
crystalline phases. This difference is often used as an estimate of 
the complexity. Already in the case of soft potentials, this gives 
rise to a number of interpretation problems; see e.g. (Binder and 
Kob, 2005; Cavagna, 2009) for a more detailed discussion. In the 
particular case of hard spheres, the situation is even worse since 
energy is irrelevant and the liquid-crystal transition is completely 
driven by entropy; hence the crystal entropy is bigger than the 
liquid one above the melting density. It follows that the con- 
figurational entropy defined as the difference between liquid and 
crystal entropy has nothing to do with the complexity; in partic- 
ular it will be negative at all densities above the melting density. 
In this paper, we will always use the word "complexity" in order 
to avoid confusion. 

12 Recall that in Fig. 1 the pressure has been rescaled. 
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FIG. 4 Schematic mean-field phase diagram of hard spheres in R , see text for a detailed description. (Left) (P, ip) diagram: 
The full black line represents the equilibrium phase diagram with the liquid-solid transition (see figure 3). The metastable 
liquid is made by a single state below ipd, while above this density it is the superposition of many glassy states. If the system 
is stuck in one of these states and compressed, it follows one of the glass branches of the phase diagram. At tpx the system 
reaches the most dense states, and if further compressed enters the ideal glass state. The pressure of the latter diverges at 
<PGCP- In the inset, the complexity, i.e. the logarithm of the number of glassy states, is plotted as function of the jamming 
density ipj. The boxes show a schematic picture of the (diV-dimensional) phase space of the system: black configurations are 
allowed, white ones are forbidden by the hard-core constraint. In the supercooled liquid phase the allowed configurations form 
a connected domain; however, on approaching ifid the connections between different metastable regions become smaller and 
smaller. Above <p k, they disappear in the thermodynamic limit and glassy states are well defined. (Right) For comparison, 
the standard glassy phase diagram for a soft potential (e.g. Lennard- Jones), specific volume vs temperature at fixed pressure, 
is reported (Cavagna, 2009; Debenedetti and Stillinger, 2001; Ediger et al, 1996). The similarity is evident if one identifies 
v — > 1/ip and T — > 1/P (i.e. one should reflect the left diagram on a diagonal line joining the upper left corner and the lower 
right corner). 



large, e.g. high-dimensional spheres or binary mix- 
tures, as discussed above. In particular, in the limit 
d —> oo we believe that the theory should be exact. 

2. In a finite dimensional system, pure equilibrium 
states cannot exist in exponential number. In fact, 
the states corresponding to ipj < p>acp can be sta- 
ble only on a finite length scale and for a finite time 
in finite dimensional systems as discussed originally 
in (Kirkpatrick and Thirumalai, 1987; Kirkpatrick 
and Wolynes, 1987a), in more detail in (Kirkpatrick 
et al., 1989; Kirkpatrick and Wolynes, 1987b), and 
more recently in (Mezard and Parisi, 2000). There- 
fore the notion of complexity makes sense only on 
a finite time and length scale. The problem of 
evaluating this length and time scales has been re- 
cently reformulated in term of a nucleation prob- 
lem (Bouchaud and Biroli, 2004; Xia and Wolynes, 
2001a, b) and is probably the most active subject 
of research in the glass transition community. The 
subject is difficult and we cannot discuss it in detail 
here. The interested reader should look to the origi- 
nal literature (Bouchaud and Biroli, 2004; Cavagna 
et al., 2007; Franz, 2005; Montanari and Semerjian, 
2006; Xia and Wolynes, 2001a,b). In the rest of this 
paper we will neglect this difficulty and assume that 
mean field states are stable on every length scale; 
however we will try to give a more precise defini- 



tion of states in Appendix A, where we explicitly 
show the origin of the difficulties in finite dimension 
and explain why our theory might be a reasonable 
approximation of the real situation. 

3. The structure of the states close to ipth might be 
more complicated than described here. In mean 
field models processes such as state crossing, tem- 
perature (density) chaos, birth and death of states, 
etc. are known to happen. This complicates the 
analysis; unfortunately at present we cannot say 
much more on this important issue. More insight 
should come from the study of mean-field hard 
sphere models such as the ones discussed in (Biroli 
and Mezard, 2001; Krzakala et al, 2008; Mari et at, 
2009). 

4. For finite dimensional systems, the transition to a 
metastable glass state is smeared and becomes a 
crossover, as seen in figure 1, see also (Ediger et al., 
1996; Moller et al., 2006). Only the ideal glass tran- 
sition has a chance to survive as a real phase transi- 
tion. For this reason the free volume fit of (Kamien 
and Liu, 2007), that describes the metastable liq- 
uid branch without assuming a phase transition, is 
not incompatible with the point of view presented 
here. 

5. When looking to actual configurations, it might 
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Density 


Definition 


Type 


Value in d = 3 


fMRJ 


Maximally random jammed configurations (Torquato et al, 2000) 


Geometrical 


~ 0.64 


<pd 


The liquid state splits in an exponential number of states 


MF, Static 


~ 0.58 


<PK 


Ideal glass phase transition - jump in compressibility 


Static 


~ 0.62 


<fith 


Divergence of the pressure of the less dense states 


MF, Static 


~ 0.64 


facp 


Divergence of the pressure of the ideal glass 


Static 


~ 0.68 


•£MCT 


Mode- Coupling transition (van Megen and Underwood, 1993) 


MF, Dynamic 


~ 0.58 




Glass transition density; depends on the compression rate 


Dynamic 


0.58-^0.62 




Divergence of the equilibrium relaxation time 


Dynamic 


~ 0.62 


<p.j 


J-point: final state of the algorithm of (O'Hern et al, 2002) 


Dynamic 


~ 0.64 


'PROP 


No general agreement on the definition 


?? 


~ 0.64 



TABLE I Summary of the relevant densities defined in the text. The values given in d = 3 (for a monodisperse system) are 
indicative, more detailed values will be given in the following. The label "MF" indicates that the corresponding concept can 
be defined in mean-field theory, but lacks an unambiguous operative definition in finite dimension. 



be very difficult to distinguish between a config- 
uration representative of a pure glass state and 
one representative of a mixture of glass and crys- 
tal (Torquato et al, 2000). Thus if one really wants 
to look at single configurations of finite systems, 
the definition of amorphous states discussed above 
is not very useful, and a definition based on or- 
der metrics may be more suitable (Torquato et al., 
2000). This leads to the concept of maximally ran- 
dom jammed packings (MRJ), defined as follows: 
one chooses a function ip that measures the or- 
der in some way, with ip = 1 corresponding to 
most ordered and tp = to most disordered con- 
figurations. Then one defines (fiMRj as the den- 
sity of the jammed configurations that minimize ip, 
see (Torquato et al., 2000) for details. This defini- 
tion may be suitable when one studies single config- 
urations, and it has been shown numerically that a 
value ifiMRj ~ 0.64 is obtained for many different 
order metrics. It is therefore important to stress 
that the phase diagram above refers to thermody- 
namic states and not to configurations. Note that 
it is likely that the local order in the configura- 
tions typical of the infinite pressure states at den- 
sity ipj depends on ipj, i.e. (ip) (tfj) is not constant; 
this does not necessarily mean that the packings 
at higher ipj cannot be considered as amorphous, 
since (ip) (ip) increases with density also in the liq- 
uid state, which is definitely an amorphous state at 
all densities. 



Before turning to the technical part, it is useful, having 
in mind this phase diagram, to come back to the behavior 
of algorithms that are used to construct jammed amor- 
phous packings, and to review some basic features that 
are observed in the structure of such packings. 



C. On the protocol dependence of the random close 
packing density 

Up to now we discussed a specific algorithm: a molecu- 
lar dynamics simulation during which the system is com- 
pressed at a given rate, the Lubachevsky-Stillinger algo- 
rithm. In a more general setting, we can consider an al- 
gorithm or experimental protocol attempting to produce 
jammed amorphous configurations of an hard sphere sys- 
tem. The algorithm stops when the system is jammed: 
the final density is a random variable depending on the 
initial data and possibly on some randomness built in 
the algorithm itself. To investigate the probability V(<Pj) 
of reaching a final density ipj, let us make the assump- 
tion that when the system is jammed, each glassy state 
contains only one configuration 13 . Then the number 
N{ifj) of jammed states with density ipj is given by 
N{ifj) ~ exp[7VE(t/?j)] where £(</?j) is the complexity 
introduced above. 

Let us define also the probability V a (<fij) that the algo- 
rithm finds one particular configuration with density ipj . 
This probability is related to the "basin of attraction" of 
configurations, i.e. to the number of trajectories of the 
algorithm that lead to the final configuration. Again, 
on the basis of concentration arguments, we expect the 
scaling V a (<fij) ~ exp[Ns a (ipj)]. 

Note that M(<pj) depends only on geometrical prop- 
erties of the configuration space of hard spheres, while 
V a {<Pj) encodes the properties of the algorithm (Krza- 
kala and Kurchan, 2007; Xu et al., 2005). The resulting 
P(<Pj) = tf(<Pj)V a (<Pj) ~ exp[AT(£(<^) + s a (<Pj))] will be 
strongly peaked around a given value ip", the maximum 
of £(<£j) + s a (<pj) (Krzakala and Kurchan, 2007). The 
important point is that (p® will depend on the partic- 
ular algorithm through the function s a (tpj), giving rise 



If this assumption is false one should change the definition of 
complexity but the following argument would remain true. See 
Appendix A for a more detailed discussion. 
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to the protocol dependence observed in experiments and 
numerical simulations (Chaudhuri et al., 2009; Xu et at, 
2005). 

Therefore it is highly non trivial to associate the phase 
diagram of Fig. 1 to the behavior of a given algorithm or 
experimental protocol. The same problem has been ex- 
tensively discussed in the context of optimization prob- 
lems, where a similar phase diagram is found, without 
reaching yet many conclusive statements. We can only 
make some remarks and conjectures, some of them al- 
ready discussed in the introduction to this section: 

1. In the static computation, infinite pressure states 
are found only in the interval [fth^GCp]- One 
cannot exclude the existence of states with very 
small volume in phase space that are not seen in the 
static computation. However, due to the smallness 
of their volume, to observe them one should design 
a very specific algorithm. Thus we can assume that 
reasonable procedures, that allow for some explo- 
ration of phase space, will always find states with 
jamming density ipj e [ip t h,<PGCp\- 

2. In the case of equilibrium dynamics (e.g. a molec- 
ular dynamics simulation with infinitely slow com- 
pression), the system is expected to equilibrate at 
each density. Thus, neglecting crystallization, we 
expect the system to follow the liquid branch up to 
<Pk and then the ideal glass branch up to focp- 

3. The Mode-Coupling theory instead predicts a di- 
vergence of the relaxation time at ^mct < fx- 
However in the context of structural glasses it is 
well known that the Mode-Coupling transition is 
avoided and becomes a crossover to a noncquilib- 
rium glass state at a compression-rate dependent 
ip g . At the mean field level, the Mode-Coupling 
transition corresponds to the point ipd where the 
liquid state breaks into an exponential number of 
disconnected states. 

4. The equilibration is particularly slow above 
</?mct ~ fdy being due to activated pro- 
cesses (Bouchaud and Biroli, 2004; Franz, 2005; 
Xia and Wolynes, 2001a, b) 14 . For moderately high 
rates, the system will leave the liquid branch very 
close to ipd and will jam at tpth- Thus the procedure 
of O'Hern et al. (O'Hern et at, 2002), that corre- 
spond to a very fast compression, should produce 
packings around tp t h- At the mean field level one 
would have ipj ~ tpth- However one should keep in 
mind (see Appendix A) that states close to ipth are 
particularly unstable in finite dimension and there- 
fore ipth is an ill-defined concept in this case. 



Point 


Densities 


Replica parameter 


Dynamic transition 


f>d, fMCT 


m = 1 


Static transition 


<fiK, <fiO 


m = 1 


Glass close packing 


<PGCP 


m = 


Threshold states 


fth, ¥>./(?)> ¥>Mflj(?) 


m = 



14 In the constant-pressure ensemble the times are proportional to 
cxp(PAV), where AV is the volume barrier to go from one state 
to another state 



TABLE II Special points of the phase diagram obtained 
within the replica computation, see figure 4. For each special 
point, we list the densities defined in table I associated to it. 
Note that (p g and <prcp depend on the protocol used, hence 
they are not associated to any special point in the replica 
computation. The identification of <p j and <pmrj with (p t h is 
tentative and probably not stricly true; still we expect these 
points to be quite close (see text for details). 



5. Remarkably, the value <Pgcp defined as the point 
where S(<Pj) vanishes is a property of the system 
that does not depend on the algorithm (at least if 
one accepts the idealizations discussed in the pre- 
vious subsections, i.e. if one neglects the existence 
of the crystal). 

6. The behavior of more complex algorithms should 
be discussed on a case-by-case basis. Thus, if one 
defines the random close packing density as the fi- 
nal density for a given protocol, its value can be 
everywhere in [pth , fee p] (and maybe outside this 
interval, for very particular protocols). Depending 
on the personal taste, one would like to identify 
fRCP with (flccp, if one is thinking to slow quasi- 
equilibrium compressions, or with tp t h ~ <Pj, if one 
is thinking to fast quenches. 

In table I we summarize the different densities defined up 
to now. 



D. Structural properties of amorphous states 

Despite the difficulty in determining the final jamming 
density for a given experimental protocol, it turns out 
that many structural properties of the final state, man- 
ifested for instance in the pair correlation function g(r), 
are roughly independent of the protocol. This suggests 
that at least some of these properties are common to 
all amorphous packings in the range ipj € \pth > Vgcp] > 
and are therefore characteristic of "random close packed" 
structures. We will be able to compute the correlation 
function g(r) for these states; therefore it is interesting to 
review the main features that are observed in numerical 
simulations. 

Detailed numerical and experimental results for g(r) at 
jamming have been reported in (Aste et at, 2005; Clarke 
and Jonsson, 1993; Donev et al., 2005a, b; O'Hern et al., 
2002, 2003; Silbert et al, 2006; Skoge et at, 2006). De- 
pending on the procedure one approaches the jamming 
point ipj from below, 6ip = ipj — ip > (Donev et at, 
2005a, b; Skoge et at, 2006) for hard particles, or from 
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above, Sip — ip — ipj > 0, for soft particles (O'Hern et al, 
2002, 2003; Silbert et al, 2006). The main features that 
are observed are: 

1. A delta peak close to r — D, due to particles in 
contact. One has g(D) <~ Sip^ 1 and the width of 
the peak oc Sip. The delta peak has the scaling form 
g(r)/g(D) = /(A), where A oc ^Sip' 1 , if <p ipj 
for hard particles (Donev et al., 2005a), with f(t) 
a scaling function independent on Sip. The area of 
this peak gives the average number of particles in 
contact with a reference one, which is found to be 
z = 2d once rattlers (particles without contacts) 
are removed; this property is called isostaticity. 

2. A square root singularity, g(r) ~ (r — D)~ a with 
a = 0.5, close to r — D. This singularity is in- 
tegrable and does not contribute to the number of 
contacts. The value of the exponent is debated: 
some claim that it is equal to 0.5 irrespective of 
the value of ipj, the procedure, etc. (Silbert et al, 
2006). However some dependence on the procedure 
has been claimed in (Donev et al, 2005a), where a 
value a <~ 0.4 has been found. 

3. A dip around r/D ~ 1.2 with respect to the liquid 
is observed. This is due to the particles of the first 
shell which are pushed toward contact (r = D) for 
Sip -> 0. 

4. A split second peak in r/D = \/3 (Donev et al, 
2005a; Silbert et al, 2006). It is not clear what is 
the exact behavior of g(r) close to this peak (Silbert 
et al, 2006), in particular if g{r) is divergent or has 

only divergent slope for r/D — > V3 . Even if g(r) 
is not divergent for r/D — > \/3, it has at least a 
jump in \/3. 

5. A similar behavior is observed in r/D = 2. 
Both features have been interpreted in (Clarke and 
Jonsson, 1993; Silbert et al, 2006) as coming from 
the network of contacts. In fact, r/D = 2 is the 
maximum distance at which two particles sharing 
one neighbor can be found, while r/D = \/3 is 
the maximum distance for two particles sharing two 
neighbors. 

6. Long range correlations (h(r) <~ 1/r 4 and S(k) <~ 
|fc| for large r and small k) are found for hard par- 
ticles (Donev et al, 2005b) and seem to be present 
also in the soft particles case (Silbert et al., 2006). 
This implies also that S(0) = 0, i.e. the packings 
are incompressible. 

7. A particularly intriguing property of jammed amor- 
phous configurations is the presence of an excess 
of soft modes, i.e. vibrational modes with very 
small frequency. This has been shown numeri- 
cally in (O'Hern et al, 2003; Silbert et al, 2005). 
In (Wyart, 2005; Wyart et al, 2005a,b) it has been 



argued that this excess of soft modes is related 
to the isostaticity property of the network of con- 
tact; moreover, a diverging length scale has been 
associated to these modes. It has also been pro- 
posed that the square- root singularity of g(r) is 
related to these modes (Wyart, 2005). This set 
of results seems to suggest a "critical" nature of 
jammed amorphous packings that is currently re- 
ceiving a lot of interest in the community (Hatano 
et al, 2007; Majmudar et al, 2007; Olsson and Tei- 
tel, 2007; Zeravcic et al, 2008). 

One of the main aims of the following discussion is to 
show that the class of packings obtained as infinite pres- 
sure glassy states share at least some of these features, 
that are common to all disordered jammed packings pro- 
duced with different protocols. We will be able to com- 
pute partially the g(r) of these states and show that it 
is consistent with the properties described above. This 
fact supports the main assumption of this paper, that the 
final states reached by typical algorithms belong to the 
class of infinite pressure glassy states. 



III. THE METHOD 

Assuming the phase diagram in fig. 4, and neglecting 
the ambiguities associated to the existence of the crystal, 
the properties of the glass phase can be computed using 
a replica method inspired by mean field models. The 
method has been described in great detail in a number 
of papers, see in particular (Mezard, 1999; Mezard and 
Parisi, 1996, 1999a, 2000; Monasson, 1995); therefore we 
will briefly sketch it here, but the reader who is inter- 
ested in the details and has no previous knowledge of the 
replica method should refer to the original papers for a 
more complete presentation. 

An alternative route to compute the properties of 
glassy states is to use density functional theory (Chaud- 
huri et al, 2005; Dasgupta and Vails, 1999; Kim and 
Munakata, 2003; Kirkpatrick and Wolynes, 1987a; Singh 
et al, 1985; Yoshidome et al, 2007). In principle the 
two methods should be equivalent (see (Castellani and 
Cavagna, 2005) for a very pedagogical discussion in the 
case of mean field models), still it seems that the replica 
method gives more accurate quantitative results and for 
this reason we will focus on this method in the following. 



A. The replica method 

The phase diagram in figure 4 is characterized by the 
existence of many glassy states at densities above ipd- At 
constant density, states are characterized by their vibra- 
tional (or internal) entropy s, defined as the entropy of 
the system constrained to be in this state without re- 
laxing toward different states. The derivative of the in- 
ternal entropy with respect to density is the pressure of 
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slope -1 



FIG. 5 A schematic representation of E(s, ip) at fixed (p. The 
behavior at small s depends strongly on the model and on the 
density (see (Krzakala et al., 2007) for a more detailed dis- 
cussion). On increasing s, Y,(s,ip) decreases and ultimately 
vanishes at a value s max (ip). The value s* is defined by 
dE/ds = — 1. For ip < <pjc, s* < s max , while for ip > ipn 
there is no solution and s* = s mox . 



the states, that is plotted schematically in figure 4. Tak- 
ing a constant ip slice of the phase diagram in figure 4, 
one will meet different states, depending on the pres- 
sure (or equivalently on the entropy). The number of 
states of entropy s at a given density ip is by definition 
Af(s) =expNT,(s,ip). 



where in the last line we performed a saddle-point 
approximation of the integral; s*(tp) is the point 
where S(s, ip) + s assumes its maximum in the inter- 
val [s m in , s ma x] ■ Indeed more compact structures have 
higher vibrational entropy (free volume) and lower com- 
plexity (their number), and the partition function is dom- 
inated by the best compromise, as expressed by the sad- 
dle point evaluation of (3). For densities if 4 < <-P < 'Pk, 
the saddle point s* falls inside the interval [s m in, s max \. 
The ideal glass transition is met at the density ipK such 
that s*{ip K ) = s max {ip K ), or equivalently E(s*,ip K ) = 0. 
Above this density, the integral (3) is always dominated 
by the upper limit of integration, therefore we have 



S(ip) 



E(s* 3 ip) + s* 

°max 

(<p) 



V < <PK , 
<p>(p K , 



(4) 



It is easy to see that (if S(s, ip) is a smooth function) 
at tpK a phase transition happens, and is manifested by 
a discontinuity in the second derivative of S(tp). Un- 
der our assumptions, close to s max (ip) we have S(s, ip) = 



£iO)(s 



■lax 

(ip)) H . ip K 

is defined by £1 = — 1 and close to ipx we have 

s*(ip) ~ s max (ip) - (1 + Ei(ip))/E 2 {cp). Thejierivative 
of S(ip) for ip - 

-T,i(<p)s' maa .(<p) 
derivative for tp 
tinuous at ifK, while it is easy to show that the com- 
pressibility jumps at ifK, as discussed in section II. B. 2 



ip K is given by S'(tp) = ^(s*,<p) ~ 
~ s' max (ip) which coincides with the 
-> ip\. Therefore the pressure is con- 



1. The ideal glass transition 

The complexity S(s, ip) (sketched in Fig. 5) is a con- 
cave function of s; it is reasonable to assume S(s, ip) to be 
a decreasing function of s, because at fixed density states 
of higher entropy correspond to more compact structures 
(in order to have more free volume) and should be more 
rare. Moreover S(s, ip) should continuously vanish at 
some value s max (ip) corresponding to the entropy of the 
best amorphous structures at this density 15 . The par- 
tition function of hard spheres at density ip is just the 
total number of allowed configurations at that density. 
In the thermodynamic limit, each relevant configuration 
belongs only to one state; and exp(A^s Q ) is the number 
of configurations belonging to the state a. Therefore one 
can write the partition function Z in the following way: 



a 



Ns a 



aw(<fl) 

dse N[^)+s] 

.(*0 (3) 



15 As already discussed one can construct denser structures by al- 
lowing a small amount of local order: we are assuming to be able 
to avoid this in some well-defined way, that will be discussed in 
the following. 



2. The replicated partition function and the (m,ip) phase 
diagram 

The basic idea of the replica approach (Mezard and 
Parisi, 1999a; Monasson, 1995) is to introduce in (3) a 
parameter m conjugated to the internal entropy of the 
states. One can think to m as a control parameter that 
at fixed density allows to select a given group of states. 

In practice this can be done by considering m copies of 
the original system, constrained to be in the same state 
by a small attractive coupling (the practical implemen- 
tation will be discussed in the following). The partition 
function of the replicated system is then 



Z m = e 



_ e NS{m,<p) ^ p Nms 



rfse A r [S( S ,^)+m S ] ^ £ N[J:(s* , V )+ms*] 



(5) 



where now s* (m, tp) is such that S(m, s) — ms + E(s, tp) 
is minimum. The introduction of the coupled replicas 
has exactly the effect of giving a weight m to the vi- 
brational entropy in (5). Only for integer m the quan- 
tity Z m has an explicit definition (and it can be eval- 
uated by direct numerical simulations), but if m is al- 
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no solution 
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FIG. 6 (Left) Schematic (m, tp) diagram: above the clustering line rrid(ip) a non trivial solution for the inter-replica correlation 
is found. This solution gives a positive complexity in the region enclosed between the lines rrid and m 3 , therefore in this region 
glassy states are present. The line m s (ip) is defined by the condition E(ra, ip) — and corresponds to the ideal glass state. The 
intersections of the line mj (m s ) with m = 1 and m — define ipd (fit) and ipth (<Pgcp), respectively. (Right) Phase diagram 
for d = 3, as obtained by solving the replicated HNC equations (section IV). Note that the static line has an unreasonable 
behavior at small m, and that the densities (fid and tpth seem too big if compared with the accepted fucr ~ 0.58 and ip j ~ 0.64 
values. 



lowed to assume real values 16 , the complexity can be es- 
timated from the knowledge of the function S(m, if) = 
ms*(m,tp) + T,(s*(m,(p),(p). Indeed, it is easy to show 
that 



s*(m,ip) 



dS(m, tp) 
dm 



E(m, (p) = S(s*(m, ip), if) — —m 
= S(m, ip) — ms* (m, ip) 



2 d[m 1 5(m, if)] (6) 



dm 



The function S(s, f) can be reconstructed from the para- 
metric plot of s*(m, ip) and S(m, ip). 

For each density, one can define a point m s (tp) as the 
solution 17 of S(m, ip) — 0, see left panel of figure 6. 
On the line m s (ip), from (5) we see that S(m,f) — 
ms max {f), then 

c ( \ - ( \ S ( rn s( ( f)> ( P) /»s 
Oglassif) = S max (f) = — . (7) 

m s [if) 

This simple prescription allows to compute the entropy 



16 This must be done by analytical continuation once the partition 
function for integer m has been computed using some approxi- 
mation. In principle, the analytic continuation might not be well 
defined, but it has been verified explicitely in mean field models 
that the procedure gives the correct results; see (Franz and Tria, 
2006; Mezard, 1999; Mezard and Parisi, 1999a; Mezard et al, 
1987; Monasson, 1995; Talagrand, 2003) for a more detailed gen- 
eral discussion of the replica method. 

17 Note that the condition S(m, tp) = is equivalent to ^q^" 1 ^ = 0, 
which correspond to the usual optimization of the free energy 
with respect to m in the 1RSB computations. 



of the glass once S{m, if) is known 18 

Generically, the entropy S(m, ip) turns out to be the 
maximum of a functional of some order parameter, that 
typically represent the inter-replica correlation functions, 
and is obtained by maximizing the functional (explicit 
examples will be given in the following). For each den- 
sity, the inter-replica correlations are non-zero only above 
some value m^ (if), and in this case the replicas are in the 
same state. In this region we obtain non-trivial values of 
Ti(m, ip), thus md(f) corresponds to the minimum value 
Smini'fi) below which there are no states and E = 0. We 
will call md{<f) clustering line, because above this line the 
space of configuration is disconnected in many clusters 
corresponding to the glassy states. The two lines md{f) 
and m s (ip) in the plane (m,<f) define a phase diagram 
which is schematically reported in figure 6 and mirrors 
the (P, f) phase diagram; in fact m is conjugated to the 
entropy, which is related to the pressure; the two phase 
diagrams are related by a Legendre transform. The lines 
m s and m^ touch at some value fTAP below which there 
are no states except the liquid one. Above iprAP, states 
are found for md(f) < m < m s (ip). When the line md{f) 
reach the line m = 1, states begin to be present in the 
liquid phase: this correspond to the point ifd- When the 
static line m s crosses m = 1, the liquid ceases to exist 



18 A very important remark (Mezard and Parisi, 1999a) is that the 
entropy of the replicated liquid S(m, tp) is analytic as long as 
m < m s (tp). In fact, the introduction of m shifts the phase 
transition that happens for tp = ipx at m = 1 to higher values 
of density for m < 1 (see figure 6). Therefore S(m, tp) can be 
computed, in the whole interesting glassy region, by analytic 
continuation of the low density (replicated) liquid entropy. 
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and the ideal glass transition is met. 

3. The equation of state of metastable glassy states 

The replica method allows to compute, for a given den- 
sity, the function T,(s,(p). But to access the equation of 
state of the metastable glassy states we would like to 
follow the evolution of each state at different densities. 
This is in general a very complicated problem already 
at the mean- field level. Therefore in order to be able 
to perform the computation, we make a very strong ad- 
ditional assumption on the phase space structure of the 
model: namely, that each state is labeled uniquely by its 
maximum possible density or jamming density <pj . This 
is the maximum value of density for this given struc- 
ture, corresponding to infinite pressure, where particles 
are in contact with their neighbors. We assume that if 
one starts from the jammed structure at ipj and slowly 
decrease the density, the particles are allowed to vibrate 
slightly around the original structure but the state main- 
tains its identity until it merges with the liquid. In other 
words, we assume that there are no bifurcations of states, 
and states can disappear only at ipj . 

If one starts in a jammed state at density ipj and 
slightly decompress the system, the state acquires a fi- 
nite entropy s(<p, <pj). We can invert this relation to get 
<Pj(s, ip). If there are no bifurcations, and we denote the 
logarithm of the number of structures with jamming den- 
sity <fj as Sj (cpj ) , we have the simple relation 

£(«,¥>) =£,•(¥>,•(*,¥>)) • ( 8 ) 

This is a consequence of our assumption 19 , because if 
states cannot bifurcate or die, their number remains con- 
stant; this implies that we can label the states equiva- 
lently by their jamming density ipj or by their complexity 

Then the procedure to extract s(<p, ipj) from Eq. (6) is 
the following: we hx a value Sj and solve 

E(m, V ) = E J -, (9) 

to get m(ip,Hj). Then we have s(<p, £j) = 
s*(m(ip, T,j),ifi). As we will see below, in the limit m — > 
we find s*{m,ip) — » — oo, i.e. the pressure diverges. The 
jamming limit then corresponds to to — ¥ 0, and the rela- 
tion between ipj and Sj is simply 

ZjiVj) = lim£(m,v?j) . (10) 



Inverting this equation we obtain ipj(Y,j), that we can 
substitute in s(ip, £j) to obtain s(tp, tpj). 

To summarize, under our assumption, in the (to, <p) 
plane of figure 4, we can draw many lines, each defined by 
S(m, (p) — Sj. Each line identifies a group of states that 
share the same jamming density ipj, which is the point 
where the line crosses m = 0, and the same complexity 
Ej . For a given line, we can compute the internal entropy 
of the corresponding states as s(ip, ipj) — s*{m{ip, ipj), ip), 
and differentiating this with respect to ip we get the pres- 
sure of the states, as drawn in the right panel of figure 4. 
Clearly if we choose £j = we recover the ideal glass 
state as discussed in the previous subsection. 



B. The molecular liquid 

The m copies are assumed to be in the same state. 
This can be implemented by constraining each particle 
of a given replica to be close to a particle of each of the 
other to— 1 replicas. The liquid is made of molecules of to 
atoms, each belonging to a different replica of the original 
system, or in other words the atoms of different replicas 
stay in the same cage. The replica method allow us to 
define and compute the properties of the cages in a purely 
equilibrium framework, in spite of the fact that the cages 
have been defined originally in a dynamic framework 20 . 
The problem is then to compute the free energy of a 
molecular liquid where each molecule is made of to atoms. 
The to atoms are kept close one to each other by a small 
inter-replica coupling that is switched off at the end of the 
calculation, while each atom interacts with all the other 
atoms of the same replica via the original pair potential. 

Note that the free energy of the replicated liquid is 
assumed to be analytic in the whole region to < m s (ip), 
as the only phase transition happens at m s (ip) when the 
complexity vanishes. Thus we can use approximations 
valid in the liquid phase to compute the free energy up 
to m s (ip). This is enough, because the free energy is 
continuous on the transition line, and therefore the free 
energy of the glass can be computed by approaching the 
line m s (ip) from below, as discussed above. Note also that 
we are interested in the regime of fairly high densities 
where we expect the cages to be small and the motion of 
atoms inside a state to be quite localized. Thus, if all the 
replicas are in the same state, inter-replica correlations 
remains strong when we switch off the coupling. Our task 
is then to compute the entropy of the replicated liquid in 
a small cage regime, when atoms in a molecule are quite 



For this reason it is often called "isocomplcxity" assumption in 
the spin glass literature. It holds for the simplest spin glass, 
the spherical p-spin glass with a single value of p (Cugliandolo 
and Kurchan, 1993), but for more general models (e.g. spherical 
spin glasses with interactions involving sets of p and q spins with 
p 7^ q, or Ising spin glasses) its exact validity is debated (Bar- 
rat et ai, 1997; Krzakala et al, 2008; Montanari and Ricci- 
Tcrscnghi, 2004); still there is a general agreement that isocom- 
plexity is approximately correct. 



Note that it has been pointed out that, already at the MCT 
level, "cages" are in fact extended objects, in the sense that 
single atoms can always hop out of their local cage, and only 
groups of many atoms can be really blocked: strictly speaking 
only when the number of atoms goes to infinity the group is 
blocked. Depending on the approximation, we will be able or 
not to take into account this effect. 
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FIG. 7 (Left) A molecule of the replicated liquid: each full (black) sphere of the original liquid is replicated m times (dashed red 
spheres), and the m copies vibrate around the reference one. (Right) Replicated entropy as a function of the order parameter 
A/Ao, where Ao is some reference value and m < 1. The full line is the mean field curve, while the dashed line take into 
account the finite dimensional nature of the system (see Appendix A) . 



close to each other. 



1. The partition function 

We start from the grancanonical partition function of 
the replicated system of molecules of coordinates x = 
(x\,-" ,%) in a volume V, and we put an harmonic 
attraction between particles in a molecule. Particles be- 
longing to the same replica interact via the hard sphere 
potential. The replicated partition function is 



oo „ 

Z m (e) = J2 zN 

N=0 JV 



d N x\ ■ ■ ■ d N x % 



■ n n xfc™ ~ x ^ 



i<j a 



i \ a<b / 



E 

N=0 



Y[z(xi)Y[x(xi,Xj) 



*<3 



(11) 



where x(x-y) = 9(\x-y\-D), x(x,y) = U. a x(x a -y a ), 
and 



z(x) = zexp — - y~](x a - x b ) 2 

\ a<b 



(12) 



in two replicas a ^ b: 



A ~ 2Nd \12( Xai ~ Xbi y 



m(m — l)d 



Xb) 



(13) 



a<b 

1 dS(m,ip;e) 
(m — l)d de 

Note that this cage radius could also be measured as the 
large time limit of the mean square displacement of a 
single system (Angelani and Fofh, 2007). We are inter- 
ested in the limit of zero coupling between the replicas. 
If there is only one thermodynamic state, the liquid, then 
the replicas will decorrelate for e — > and A(0) = oo. On 
the contrary, if the are many stable states, an infinites- 
imal coupling will be enough to send the replicas into 
the same state; this will produce a correlation between 
the replicas and A will be of the order of the cage radius 
inside one state. This phenomenon can be better un- 
derstood in term of the Legendre transform of S(m, (p; e) 
with respect to e, which is a function of A defined by 21 



S(m, ip; A) — min [<S(m, ip; e) + (m — l)cL4e] 



(14) 



and is schematically drawn in the right panel of fig- 
ure 7. As (to — l)de = dS (™^' A } ; the stationary 
points of S(m, ip; A) correspond to zero coupling and in 
fact S(m,(p) = maxyi ^(m, ip; A). In the liquid phase, 
S(m, p; A) will have a single maximum in A — oo, while 



It is clear that the derivative of <S(to, ip; e) 

Y 1 In Z m (e) with respect to e gives the average cage 
radius, which is defined as the average distance of atoms 



Due to a global m — 1 factor, in the following discussion minima 
and maxima are exchanged for m < 1. 
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when glassy states are present, S(m, ip; A) will have a 
secondary maximum at finite A. Sec (Mczard and Parisi, 
2000) for a much more detailed discussion of the Legen- 
dre transform with respect to e. 

We will discuss in the following some approximation 
schemes to compute Z m (e). For the moment, let us try 
to give a "pictorial" representation of this method. A 
way to visualize the partition function Z m (e), Eq. (11), 
is that the original system is represented by the reference 
coordinates xu of the first replica; around each reference 
particle, m—1 other particles vibrate on a small scale ~ 
\j\ft- They are represented in Fig. 7 as dashed spheres. 
The interaction is such that dashed spheres belonging to 
the same replica cannot overlap. 

Alternatively, the system can be thought as a "molecu- 
lar liquid" in which each molecule is built by the m repli- 
cated particles vibrating around their center of mass X = 
mT 1 J2 a x a . The entropy of the replicated system is given 
in Eq. (6) by S(m,ip) — Y,(s*(m,(p),ip) + ms*(m,ip). 
In this picture, the entropy of the molecular system, 
S(m, ip), is given by the entropy of the centers of mass, 
E(m, <p), plus the contribution due to the vibrations of 
the m particles in a molecule around the center of mass, 
which is m times the vibrational entropy s*(m, <p), i.e. 
the volume of a typical cage of the centers of mass sys- 
tem. Note that for generic m, s* slightly depends on m 
since the m replicas are not independent. 

This representation of Z m (e) shows that a consistent 
way to remove the crystal state in this computation is 
to describe the system of the centers of mass by a low 
density virial expansion, as it will be done in the following 
(see Appendix B). In other words we are assuming that 
the centers of mass X represent structures that are typical 
of the liquid state. This is our operational definition of 
amorphous states and for this reason in our computation 
partially crystallized states do not appear. 



2. Correlation functions 

The replicated liquid is characterized by the density 
pa{x) — (J2i ^i x ia — x)) of each replica and by the cor- 
relation function 

Pab(x, y) = (^2S(x ai - x)S(x bj - y)^j . (15) 

where the sum is over all ij is a ^ b and over i ^ j 
if a = b. If there is symmetry between the replicas, 
p a {x) = p and the two relevant correlations are the intra- 
replica correlation g(x,y) = p aa (x,y)/ p 2 and the inter- 
replica correlation g(x, y) = p a ^b{x, y)/p 2 - In the follow- 
ing we will show that g(x,y) has a quite different shape 
in the liquid and in the glass at high pressure. Hence, 
we will denote by ga{x,y) the pair distribution function 
in the glass phase, while we keep the notation g(x,y) 
for the one of the liquid. To avoid confusion, we stress 
that g{x,y) and gc(x,y) refer to the same observable in 



different phases, while g(x,y) and g(x,y) are different 
observables. 

If there are no correlations between different replicas, 
g = 1. This happens for m < ma (ip), where no glassy 
states are present. In the region m > rrid (<p), there 
are glassy states. Each state a is characterized by a 
frozen density profile p a (x) and by its correlation func- 
tion p a (x,y). The density profile can be computed, in 
principle, as the minimum of a suitable density func- 
tional F[p(x)] (Chaudhuri et at, 2005; Dasgupta and 
Vails, 1999; Kim and Munakata, 2003; Kirkpatrick and 
Wolynes, 1987a; Singh et at, 1985; Yoshidome et at, 
2007). 

We are interested in a verage s over the states. Due to 
translational invariance, p a (x) = p, and for the intra- 
replica correlation we have simply 

g{x - y) = p~ 2 p a (x,y) , (16) 

where the over line denotes average over the states. 

We assumed that the coupling between different repli- 
cas is able to force all them to be in the same state. 
Apart from that, in the limit of vanishing coupling no 
other correlations are induced. Thus 

g(x ~y) = p~ 2 p a {x)p a {y) . (17) 

3. Nonergodicity factor 

Eq. (17) allows to relate the inter-replica correlations 
to the so-called nonergodicity factor of Mode-Coupling 
theory. We sketch here the connection but the reader 
is referred to (Bengtzclius et ai, 1984; Gotze, 1999; van 
Megen and Underwood, 1993) for more details. In Mode- 
Coupling theory it is usual to work in Fourier space; the 
Fourier-transformed density reads 

p a (q)= J dxe i ^p a (x) = (j2^ 9Xi ^ - (18) 

where (*) a denotes thermal average in state a. The main 
object of Mode-Coupling theory is the coherent normal- 
ized scattering function 

F[q > t] = nW) ^E e49Mt) ^' (0)1 ) - ( 19 ) 

where S{q) = iV" 1 ^ e^" 3 ^ = 1 + ph(q) is the 
structure factor, which is related to the Fourier trans- 
form of h(r) — g(r) — 1 (Hansen and McDonald, 1986), 
and Xi(t) is the position of particle i at time t assum- 
ing that the position at time t = has been extracted 
from the Gibbs distribution (Mode-Coupling theory is a 
theory of the equilibrium dynamics of glasses). By defi- 
nition, F(q, 0) = 1; in the glass phase F(q,t) develops a 
plateau that becomes infinite at the Mode-Coupling tran- 
sition (Pmct- Hence for ip > ipmct the long-time limit 
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FIG. 8 The three diagrams contributing to order four in den- 
sity to the replicated free energy functional. Each vertex is 
associated to a density p a , while each line is associated to a 
hab(xi, Xj) = g a b(xi, Xj) — 1 factor; a, b = 1, • • • , m are replica 
indeces. The value of the diagram is obtained by integrating 
over the Xi and summing over the replica indeces. Only the 
leftmost diagram is kept in the HNC approximation (Hansen 
and McDonald, 1986); the other two are two- lines irreducible. 



and, taking into account that the initial condition (hence 
the initial state a) has been extracted from the Gibbs 
distribution, one finally obtains 



p a (q)p a (-q) = ph(g) 

NS(q) S(q) 



(22) 



where h(q) is the Fourier transform of h(r) = g{r) — 1. 



IV. THE REPLICATED LIQUID: HNC EQUATIONS 



of F(q,t) is non zero; this is called nonergodicity factor 

u 



f q =\im F(q,t) . 



(20) 



In the replica interpretation, the fact that F(q,t) does 
not vanish in the long time limit (or in other word, that 
density fluctuations cannot relax completely) signals the 
appearance of metastable states. Even for t — > oo, the 
system cannot escape from the metastable state in which 
it was at time t = 0. However, it can decorrelate inside 
the state. Therefore, in the long time limit one has 



yygM«)-sj(Q)] 

\ ij I 

= pa{q)pa(-q) , 



iqxi(t) 



a \ J 



4^ 



The simplest way to compute the property of the repli- 
cated liquid is the following. We consider the system 
of the m replicas as a m-component mixture, which is 
then described by the number density p a of type a par- 
ticles, and by the correlation function g ab (x,y) which is 
the probability of finding a particle of type b in y given 
that there is a particle of type a in x. 



A. Replicated HNC equations 

The entropy of the replicated liquid can be expressed 
as a functional of these quantities (De Dominicis and 
Martin, 1964; Hansen and McDonald, 1986; Morita and 
Hiroike, 1961): 



S[p a ,gab(x,y)} =~\y2 f dxdyp a p b {g ab (x,y)\ng ab (x,y) - g ab {x,y) + g a b{x,y)p(j) ab {x,y) + 1] 

ab J 

^ (—1)™ 

— V p a [In p a — 1] — - Tr[p/i]™ — {two-lines irreducible diagrams} , 



(23) 



and the g a b(x,y) have to be determined by maximizing the entropy. In the expression above we assume that the 
inter-replica coupling has already been sent to zero and look for non-trivial solutions for the inter-replica correlation. 
As an example, the diagrams contributing at order four in density are given in figure 8. 

The simplest approximation that allows to obtain a treatable functional amounts to neglect two-lines irreducible 
diagrams; in this way one obtains the HNC equations for g ab '- 



ln 5a6 (x, y) + l3(f> ab (x, y) = h ab (x, y) - c ab (x, y) 
where c is defined by the OZ relation 



Kb(x,y) = c ab (x,y) J dzh ac (x,z)p c c cb (z,y) . 



(24) 



(25) 



The interaction potential is <fi a b{x,y) = 4>(x — y)S ab , where <fr(x) is the hard core potential. Using translational 
invariance, and a replica symmetric structure, we have p a = p, g aa (x,y) — g(x — y) and g a ^b(x,y) — g(x — y). The 
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replicated entropy becomes 22 



S(m,ip) 



S_ 

N 



+ 



dr {mg{r) [In g(r) — 1] + m(m — l)g(r)[ln g(r) — 1] + m 2 + m(j>(r)g(r)} — m(lnp — 1) 



2 2 

where /1(g) are the Fourier transforms of h(r), h(r). The HNC equations can be written as 

\ng(r) = -<f>(r) + W(r) , 
In g(r) = W(r) , 

where W(r), W(r) are defined by their Fourier transforms 



W(q) 
W(q) = 



1 p[h(q) + (m - l)h(q)} 2 + m-l p[h(q) - h(q)} 2 



mi + p [h(q) + (m-l)h(q)} m 1 + p[h(q) - h(q)} 
1 p[h(q) + (m - 1)%)] 2 1 p[h(q)-h(q)Y 



m 1 + p[h(q) + (m - l)/i(g)] m 1 + _ h ( q )} 

Using Eq. (6), the free energy and complexity of the states are 
dS(m) 



s*{m) 



dm 



- P - J dr{g(r)[\ng(r) - 1] + (2m - l)flf(r)[liiflf(r) - 1] + 2m + <f>(r)g(r)} - lnp + 1 

+ hS^ Hl+p{h{q) - lm+ 



ph(q) 



l + p(h(q) + (m-l)h(q)) 



_ p%)+ ^i! + ( 2 --y%) 2 }, 



S(m) = -m 29(5(m)/m) - Pm2 



9m 



| dr{g(r)[lng(r)-l] + l}-±J ^{ ln[l + - %))] 



-ln[l + p(%) + (m-l)%))] 



mph(q) 



1 + + (m - 1)%)) 



+ 



m 2 p 2 h(qY 



(26) 



(27) 



(28) 



(29) 



The advantage of this formulation is that equa- 
tions (27) are relatively easy to solve numerically and 
give direct access to both the thermodynamic property 
of the glass (entropy and pressure) , the structure factor 
g(r) and the non-ergodic parameter g(r). 



these calculations for illustrative purposes 23 . 



1. Phase diagram 

The (m, ip) phase diagram for d = 3 is reported in 
figure 6. The values of (fx = 0.63 and = 0.619 are 
reasonable, even if ipd is a bit too large if compared to 
(Pmct ~ 0.58 (van Megen and Underwood, 1993). Also 



B. Results 



These equations were used in (Cardenas et at, 1998, 
1999; Mezard and Parisi, 1996) to compute the proper- 
ties of amorphous states of hard spheres. We reproduced 

22 We can set f3 = 1 as temperature is irrelevant for hard spheres. 
Note that the term f dr<f>(r)g(r) = 0. 



The HNC equations were solved by an iterative Picard scheme, 
using an initial Gaussian guess for c(x). We used a grid defined 
by imposing a cutoff r < L = 8D and discretizing space with a 
step a = D/128. We checked the stability of the reported results 
by doubling the cutoff and inverse step. 



22 



the value of ipth ~ 0.67 that one can guess from the 
extrapolation of nid to m — is bigger than the accepted 
value ip. j <~ 0.64 from numerical simulations. 

Moreover, one immediately notices that the static line 
does not seem to extrapolate to zero at reasonable densi- 
ties 24 . This was already observed in (Mezard and Parisi, 
1996) and is one of the major problems of the HNC ap- 
proximation. Before discussing this issue in detail, let us 
discuss the results for the correlation functions. 



2. Correlation functions 

Clearly in the liquid phase m = 1 and the correlation 
function of the liquid coincides with the one obtained 
within the non-replicated HNC approximation. The re- 
sult for ga{r) along the line m s (ip) (i.e. for the ideal 
glass) is reported in figure 9. Many interesting features 
expected from numerical simulations (see section II. D) 
are observed: first of all, gc(r) develops a delta peak 
close to r = D, whose shape (inset of left panel) is quite 
similar to the one observed by Donev et al. (Donev et al, 
2005a). Moreover, a dip for r ~ 1.2 is developed, due to 
particles of the first shell getting closer to the reference 
one. Another interesting feature is the jump that devel- 
ops in r = 2D and the cusp in r = 3D, which are in fact 
related to the delta peak in r = D. Finally, we observe 
that on increasing the density go (r) develops long ranged 
oscillations for large r, even if a systematic study of these 
correlation would require the investigation of smaller val- 
ues of m at which the approximation gives inconsistent 
results. 

As for g(r), the most interesting feature is the peak 
close to r = that describes the particle of replica b 
which is in the same molecule of the reference particle of 
replica a. Its shape has a clear scaling form on decreas- 
ing m: the cage radius, i.e. the scale on which the delta 
peak is rounded, decreases and goes to for m — > 0, 
as anticipated in previous sections. The behavior for 
larger r is very similar to the one of g{r), "rounded" 
by the fluctuations inside each molecule of particles of 
different replicas. Finally, in figure 10 we report the 
nonergodicity factor f q , Eq. (22), which is a central ob- 
ject of Mode-Coupling theory. The shape of f q obtained 
with replicas is qualitatively very similar to the one ob- 
tained within MCT (Bengtzelius et al, 1984; Gotze, 1999; 
van Megen and Underwood, 1993) and experiments (van 
Megen et at, 1991), although from the quantitative point 
of view marked differences are present. 



It is worth to note at this point that on increasing the density 
long range correlations seems to develop and short range singu- 
lar behavior emerges in g(r), as expected from numerical simula- 
tions. Thus one is forced to increase the cutoff and inverse step 
used in the discretization of the HNC equations. Wc checked 
that in the range reported in figure 6 the resulting m s , are 
not affected by discretization corrections. 



C. Discussion 

The replicated HNC equation already gives interesting 
qualitative indications on the phase diagram and correla- 
tions functions, and is not too bad from the quantitative 
point of view. However, the value (fd — 0.619 is too large 
compared with numerical estimates, and the nonergodic- 
ity factor is quite far from the measured values. Another 
unsatisfactory feature of the HNC approximation is that 
it gives a complexity E of the order of 0.01, which is 
two orders of magnitude smaller compared to what is ob- 
served in simulations and experiments, where E ~ 1. As 
for the correlation of the glass, this approximation misses 
the peak at r = \/3D (which should be encoded in the 
diagrams that have been neglected) and the square-root 
singularity close to r = D, but reproduces well the other 
characteristic features of jammed disordered packings. 

Unfortunately, the results of the HNC approximation 
become definitely bad at small values of m, i.e. deep in- 
side the glass phase on approaching jamming. This can 
be understood quite easily by inspecting the diagrams 
that have been neglected in the HNC approximation, al- 
ready at the lowest non-trivial order in density; these are 
shown in figure 8. The argument goes as follows: as dis- 
cussed above, when the pressure is high (or m is small) 
the cage becomes very small, and the peak in g(r) close 
to r = approaches a delta function. Consider the con- 
tribution to the diagrams in figure 8 when the replica 
indeces are all different, a ^ b ^ c ^ d; then on each link 
there is a factor g(xi — Xj). We focus on the small X^ X j 
behavior; then g(xi — Xj) <~ 6(xt — Xj) and the leftmost 
diagram has a contribution 

J dx 1 dx2<ix' i dx i 5{xx — X2)5{x2 — x^)x 
5(x 3 — x 4 )5(x 4 — x\) ~ 5(0) 

since three delta functions are enough to constrain the 
four integration variables to be close to each other. The 
notation 5(0) here is a formal way to indicate that the 
contribution above will diverge roughly as the maximum 
of the peak of g(r). Now, the other two diagrams in fig- 
ure 8 have more links; performing the same computation 
as above, we find that they diverge as 5(0) 2 and <5(0) 3 . 
But these diagrams are not included in the HNC and this 
makes the approximation bad at high pressure. Actually, 
this reasoning shows that the HNC is keeping the less 
divergent diagrams at each order, while the maximally 
divergent diagrams are the completely connected ones, 
at least as far as the contribution with different replica 
indeces are concerned. Therefore to study the jamming 
limit we need to treat correlations inside a molecule in a 
more appropriate way. We will discuss in next section a 
different approximation scheme that fixes this problems 
by taking into account all the correlations between differ- 
ent replicas (in some sense, resumming the contribution 
of completely connected diagrams). 

The argument above is supported by the fact that we 
were not able to derive analytically the limit m — > of the 
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FIG. 9 gaif) and g[r) along the ideal glass line m s {<p) at three different densities, computed using the replicated HNC equations 
(section IV). Inset: scaled plot of gair) / <?g(-D) vs p(r — D)/D and g(r)/g(0) vs pr/D, where p = 1 + 4:Lpgc(D). 




FIG. 10 Nonergodicity factor, defined in Eq. (22), com- 
puted using the replicated HNC equations (section IV) at 
ip = 0.62 ~ (fid and tp = 0.63 ~ <£k ■ 



HNC equations, and we have a strong feeling that they 
become pathological in this limit. The same problem is 
present in the limit d — > oo, and should be connected to 
the fact that the cage radius is very small in this limit, as 
we will see below. In summary, the HNC approximation 
does not work for small cage radius: multiple correlations 
functions in the replica sector become important. 



V. THE REPLICATED LIQUID: EFFECTIVE 
POTENTIALS 

A way to compute the replicated free energy for a 
system of Lennard- Jones particles in a more accurate 
way, by taking into account multiple correlations inside 
a molecule, has been described in (Mezard and Parisi, 
1999b). The idea is to write the free energy as a function 
of the cage radius A and expand systematically in powers 



of A, which is assumed to be small. Thus we expect this 
method to work well in the dense region and for m ~ 0, 
where the HNC approximation fails. The method is suc- 
cessful but cannot be extended straightforwardly to hard 
spheres, because at some stage in (Mezard and Parisi, 
1999b) it was assumed that vibrations were harmonic, an 
approximation that clearly breaks down for hard core po- 
tentials. In this section we will discuss a general method 
that allows to map the replicated free energy onto the 
non-replicated free energy of a liquid of particles interact- 
ing via some effective potentials to be computed below. 
The method is very similar in spirit to the self-consistent 
phonon theory used in (Hall and Wolynes, 2003, 2008; 
Stoessel and Wolynes, 1984) to study the stability of a 
glass state. This will allow to derive in a simple way the 
small cage expansion for the case of hard spheres. More- 
over, the correlation function of the glass turns out to be 
the correlation function of the effective liquid, thus sim- 
plifying a lot its computation. Similar ideas have been 
recently used in (Stevenson et ai, 2008). 



A. Entropy as a functional of the single-molecule density 

We wish to compute the entropy as a functional of 
the order parameter A, the cage radius, whose expected 
behavior has been schematically illustrated in figure 7. 
Instead of expanding directly the partition function (11) 
for large e, and then Legendre transform with respect to 
e as in Eq. (14), we will first use standard liquid theory 
to Legendre transform Eq. (11) with respect to the full 
function z{x). In this way we obtain the entropy as a 
functional of the single-molecule density 



P{x) = ( ^S(x - Xi) 



(31) 
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where 8{x — Xi) 
function x(x, y): 

S[p(x),x(x, 



Ila ^( x a ~ x ai)i and of the interaction 



dxp(x)[l — lnp(x)] 



^ [a class of diagrams] . 



(32) 



The class of diagrams contributing to (32) is defined pre- 
cisely in (De Dominicis and Martin, 1964; Hansen and 
McDonald, 1986; Morita and Hiroike, 1961) but is not 
important here. What is important is that a diagram T> 
represents an integral of the form 

v =\ /rip(^n(xw-i), (33) 

where i are the vertices of the diagram, I = (i < j) arc the 
links, and to lighten the notation we defined p(i) — p(xi) 
and x(^) = x{xi,Xj). S is the symmetry factor of the 
diagram, i.e. the number of equivalent relabelings of the 
vertices of the diagram. 

The reason why we start from (32) instead of (11) is 
that the single-molecule density p(x) is directly related to 
the order parameter A we want to study. In particular, 
we can make a simple Gaussian ansatz 25 



p(x) 



pm' 



-d/2 



(2TTA)( m - 1 ) d / 2 



e 2 



(2irA) md / 2 

= PJ dX U^ A)d/ 2 

= pj dX]\ lA {x a -X) 



dXe-^^- {Xa - x)2 
hi^-x) 2 



(34) 



where ja{x) is a normalized Gaussian with variance A; 
we have then J dxp(x) = pV = N, therefore p has the 
interpretation of the number density of molecules (note 
that having Legendre transformed we are now working 
at fixed N). 

The Gaussian ansatz for p(x) is not equivalent to the 
Gaussian ansatz (12) for the single molecule activity; in 
fact p{x) is a sum of diagrams containing z(x) and the 
interaction function, so it is Gaussian at lowest order in 
e but has corrections coming from higher orders. The 
corrections are singular because of the singularity of the 



25 This simple ansatz assumes that all particles have the same cage 
radius \f~A. However, at very high pressure, it is well known 
that most of the particles arc immobile, while a small fraction 
(typically ~ 5%) can vibrate in a cage that remains finite even 
at infinite pressure. These particles are called "rattlers" in the 
literature on jamming. Our ansatz completely neglects rattlers. 
One could think to improve it by introducing two cage radii, one 
associated to jammed particles and the other to rattlers: we did 
not explore this possibility. 



interaction. For this reason it is convenient first to Leg- 
endre transform with respect to z(x) and then to make 
the Gaussian ansatz for p(x): the resulting small A ex- 
pansion has a better behavior. 

The Gaussian form (34) allows to compute exactly the 
ideal gas term of the entropy: 



J dxp(x)[l-lnp(x)] = 



1 - In p - ^ (1 - to) ln(27rA) - ^ (1 - to - In to) 
1 -lnp + S h 

arm 



(35) 



where we defined 



S h arm(m,A) = ^ (m- 1) ln(27r A)+^ (m-l+ln to) . (36) 

We want to find a simple way to rewrite the replicated 
entropy exploiting the fact that A is small. 



B. The effective non-replicated liquid 

Before proceeding to the formal computation, it is bet- 
ter to understand intuitively what will be the result. To 
this aim, note that we assumed that the vibrations of 
the to copies of the particles, x — (xi, • • • , x m ), are de- 
scribed by the Gaussian distribution (34), corresponding 
to harmonic vibrations. 

The width A is a variational parameter and we will 
maximize the entropy with respect to it at the end; for 
the moment we assume that A is small. The idea is to 
choose a replica (say replica 1) as a reference and con- 
sider the vibrations of the other to — 1 particles around 
the reference one. At the zero-th order in A, the to — 1 
copies essentially coincide with the reference one, and 
iS(m, ip\ A) = iV _1 5[p(x), x(x, y)] is given by the en- 
tropy S(ip) of the non-replicated liquid (corresponding 
to replica 1), plus the free energy of to — 1 harmonic 
oscillators of spring constant A: 

S<°> (to, A) = S{<p) + S harm (m, A) . (37) 

We see indeed that the ideal gas term (35) corresponds 
to Eq. (37) where S(ip) has been approximated by the 
ideal gas contribution. 

A first order approximation is obtained by considering 
the effective two-body interaction induced on the parti- 
cles of replica 1 by the coupling to the to — 1 copies. Wc 
consider two particle x\ and y\ of replica 1 belonging to 
two molecules x and y. Each particle of a given replica a 
interacts with the other particles of the same replica via 
the hard core potential <j){r). The effective interaction 
between particles in replica 1 is obtained by averaging 
this interaction over the probability distribution p(x) of 
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the two molecules x and y: 



r. m 

■>"&-v^= J dx 2 , m dy 2 , m p- 2 p(x)p(y)l[e-^-y« 



= e -<t>(xi-yi) I TT e -<t>(x a -y a ) 



n 

\o=2 



xi, yi 



(38) 



As an example, the potential <p e ff(r) for hard spheres in 
d = 3 is reported in Fig. 11 for some values of A and m 
(see Appendix C for its calculation). A first order ap- 
proximation to 5(m, (p; A) is then obtained by substitut- 
ing to the entropy S(ip) of the simple hard sphere liquid 
the free energy of a liquid of particles interacting via the 
potential <peff(r)- 

SW(m,ip;A) = -(3F[<p;<t> eff {r)] + S harm {m,A) . (39) 

It is now evident that we can obtain better approxima- 
tions of the true function S(m, ip; A) by considering also 
the three body interactions induced on particles of the 
replica 1, and so on. 

This procedure can be justified more formally by in- 
troducing a diagrammatic expansion in powers of A and 
resuming a class of diagrams: this is discussed in detail in 
Appendix B. The result is the following: the replicated 
entropy (32) can be exactly rewritten as 

S(m,<p;A) = N-'Sipix),^^)} 

(2) ('3') 

= S harm (m,A) - 0F eff [<p;(t> e ff,(t>y f ,(t>l/ f ,---] , 

(40) 

where the effective potentials 4>^j depend on both A 
and m. The correlation function of the glassy states as 
function of A and m is just the correlation function of 
the liquid described by F e j / . 



The effective potential 



is constructed as follows: 



one constructs all the possible connected diagrams with 
n links and such that one pair of vertices is not connected 
by more than one link. Then one numbers the vertices 
Xi, X2, ■ ■ ■ ■ The effective potential (jr™^ for one chosen di- 
agram depends on the variables x\, x 2 , • • • (whose num- 
ber depend on the diagram). For n = 1, we have only one 
possibility (j) e ff depends only on r — x\ — x 2 by trans- 
lation invariance and is given by (38). For n = 2 also we 
have only one possibility where the two links share one 
vertex, the resulting potential depend on r = x\ — x 2 and 
s = x\ — x 3 (denoting by x\ the common vertex) and we 
have 



,(2) 



ff (x-y,x-z) 

<nr= 



,-4>(x a -y a ) p -<t>{x a -z a )\ 



(41) 



<n: =2 



».)) 



where the brackets denote averages similar to the one in 
(38), see Appendix B for details. There it is also shown 
that the correlation function (feM of the glass is given 
by the correlation function g e ff{r) of the effective liquid. 



C. Properties of the effective potentials 

The effective potentials have the following general 
properties, that are discussed in Appendix B: 



1. If for at least one of the links the distance r = 
| X<1 ~~ Xj | IS such that \i — D\ >• >/A, the potential 
vanishes exponentially as e~( r ~ D ) / A . 



2. For the reason above, we can argue that they con- 
tribute 0(A n / 2 ) to the free energy, where n is the 
number of links in the potential. This is because 
the potential is non-vanishing only if, for all the 
links, \r — D\ <~ \/A, and this region has volume 
0{A n ' 2 ). 

3. One can show that the leading order of the po- 
tentials vanishes if all the unit vectors of the links 
are orthogonal. Then, one can argue that in the 
limit d — > oo where the links are almost always or- 
thogonal, the potentials for n > 2 give vanishing 
contributions. 

4. All the potentials are 0(m — 1) for m — > 1. Thus 
in this limit they can be treated as perturbations 
of the hard spheres liquid. 

5. In the opposite limit A — > with A — am, (jam- 
ming limit, see below) the potentials tend to be- 
come delta functions around r = D. 

In the following we will use these properties to derive the 
glass equation of state and its structure factor as function 
of the dimensionality. 

Before concluding this section it is worth to remark 
that in the case of a smooth potential it might be prefer- 
able to use as reference positions the center of masses 
instead of the coordinates of replica 1. The expansion 
of Appendix B can be carried out similarly in this case 
and leads to analogous expressions for the effective po- 
tentials. The difference is that in the case of a smooth 
potential, the small cage expansion starts with a term of 
order A and one can show that if the center of masses 
are chosen as reference positions, the n-body potential is 
of order A" 1 ^ 1 . Conversely, if one uses replica 1 as a ref- 
erence, all the potentials are of order A. The advantage 
of using replica 1 as a reference is that the correlation 
function of the effective liquid is directly the correlation 
function of the glass. In the case of hard spheres the ex- 
pansion is well behaved because of the properties above. 
One should keep in mind that depending on the problem 
at hand, one expansion might be better than the other. 



VI. THE LIMIT OF LARGE SPACE DIMENSION 

We will begin by the study of the limit d — > oo, because 
in this limit all the expressions simplify a lot and we are 
able to obtain a consistent solution in all the regions of 
the phase diagram. Moreover, this limit is interesting 
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FIG. 11 The non-hard core part of the effective potential 
4>eff{ r ) — <M r ) — — Ml + Q( r )]> Eq (38), as a function of (r — 
D)/Vtl for r > D for different values of m = 0.01, 0.1, 0.5, 3 
(from bottom to top) in d — 3. Note that the range of the 
potential is 0(y/A) and that it is attractive for m < 1 and 
repulsive for m > 1. The strength of the potential diverges 
for m — > 0. 

because metastability effects should be less important in 
high dimension: the limit d — > oo is a kind of "mean 
field" limit where the surface and the volume are of the 
same order of magnitude and nucleation becomes almost 
impossible. As we briefly discussed in the introduction, 
the limit d — ► oo has also interesting applications in the 
digitalization of signals, see (Conway and Sloane, 1993). 



so that the above expression is well defined. If 2 p is not 
exponentially large, the non-ring diagrams are conjec- 
tured to give only exponentially small corrections (Frisch 
and Percus, 1999), and one can show that the last term 
is exponentially small (Frisch and Percus, 1999). There- 
fore, if 2 d ip does not grow exponentially with d, S(ip) is 
given by the ideal gas term plus the first virial correction 
(i.e. by the Van der Waals equation) up to exponentially 
small corrections: 

S(ip) = l-]np-2 d - 1 ip + 0(e- d ) , 
g(r)= X (r)[l + 0(e- d )}. 

In (Parisi and Slanina, 2000) simple equations for the pair 
correlation function g(r) were introduced, and solved in 
the limit d — > oo; it was shown that one obtains the 
same result, up to <Pmax> meaning that Eq. (44) might be 
a reasonable analytic continuation of the liquid equation 
of state up to p^ax which is exponentially larger than 
2 d . At this value of the density a pole develops at finite 
q that seems to correspond to a liquid instability (the 
Kirkwood instability (Frisch and Percus, 1999)). 

The physical meaning of this instability is not clear. 
Although very interesting from the mathematical point 
of view (maybe also in relation to the problem of find- 
ing the most dense lattices (Parisi, 2008)), this question is 
not physically relevant in this context. We will show later 
that the glass transition indeed preempts this instability 
that is therefore in a non-physical region of the density: 
this system becomes unstable toward replica symmetry 
breaking at a density <p ~ 2~ d d, before reaching the Kirk- 
wood instability. 



A. The liquid in d — > oo 



B. The effective liquid: the Baxter model in d — > oo 



The problem of computing the entropy of the hard 
sphere liquid for d — > oo was addressed in (Frisch and 
Percus, 1999; Parisi and Slanina, 2000), where the same 
result was obtained in two independent ways. In (Frisch 
and Percus, 1999) it was shown that the ring diagrams 
dominate the virial series order by order in p for large d. 
The resummation of these diagrams gives 



%>) = l-lnp+£ J dff(r) 
dq 



2p J {2Tr) d 



]n[l-pf(q)]+pf(q) + 



(p/(?)) S 



(42) 



where f{q) is the Fourier transform of the Mayer function 
f(r) = x(r) - 1, 

f(9) = -V d r(~ + l) " J<(q) . (43) 

The key observation is that, up to a density (p^ax = 
exp[d(l — ln2)/2]/2 d , the logarithm does not have poles 



To compute the free energy of the replicated liquid, we 
will neglect all the potentials but the two-body one. This 
is because as discussed above we have indications that 
the corrections coming from the many-body potentials 
vanish in this limit. However this point deserves a more 
careful investigation. 



1. Two-body potential in large dimension 

The calculation of the two-body potential was already 
done at first order in (Parisi and Zamponi, 2005); in Ap- 
pendix C it is carried out at any order. The resulting 
effective potential, Eq. (38), is reported in figure 11 for 
d = 3. Note that at fixed m the effective potential is a 
function of (r — D)/\/4A only. 

In the limit d — > oo the correlations of the liquid vanish 
and its g(r) approaches a step function 6(r — D); on the 
contrary, as we will show below, the cage radius in the 
relevant region is very small (of the order oil/d). Then, 
as the two-body potential vanishes on a scale y/A, it is 
very reasonable to approximate it with a delta function 
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FIG. 12 The function Ti(A). 



in r = D: this leads to the Baxter model (Baxter, 1968) 
where 



-*'"M= x (r)[l + Q(r)] 

-^ X (r)[l + DG m (A)V d (l)S(r-D)] , 



(45) 



where G m (A) is related to the integral of Q(r) by 
Eq. (C4) of Appendix C: 



G m (A) 



V d (D) 



drr d - l Q(r) 



(46) 

We will show below that the natural scale for A in this 
limit is A — d 2 A/D 2 which is of order 1 at the glass 
transition. With this scaling it is possible to compute 
the leading expression in the limit d —> oo, see Eq. (C32) 
in Appendix C: 



lim G(m, D Ad ) = Q m {A) 

d—*oo 



/'CO 

/ dye y 


r / -~-\ m -i 


> — OO 





(47) 



where 6(t) = |[l + erf(t)]. 



2. The Baxter liquid in large dimension 

We have then to solve the Baxter model in the limit 
d oo. To do this it is enough to observe (Parisi and 
Slanina, 2000) that the Fourier transform of 6(r — D) 
and of 9(D — r) coincide at the leading order in the limit 
d — s- oo. Using the results of (Parisi and Slanina, 2000) 
it is easy to show that 



feff(q) = U-G m (A))f( q ) 



(48) 



where / e // is the Mayer function corresponding to the 
Baxter potential (45). Note that this relation is exact for 



q = 0. It is possible to show (e.g. by direct numerical 
computation) that, for all m and A, Q m (A) < 1; therefore 
the coefficient 1 — Q is positive and the only difference 
between the Baxter liquid and the hard sphere liquid is 
in a 0(1) renormalization of the density. Therefore, up 
to a density (fi^ a a f er = ip^ x /(l-G m {A)), the free energy 
of the Baxter liquid is given by 



-f3F Bax ter(m, V ;A) = 1 - lnp - 2 d ~ V (l - Qm{A)) . 

(49) 

We will see that the replica instability happens for cp < 
Pmax er > so this expression will be enough for our pur- 
poses. The replicated entropy is finally given by substi- 
tuting (49) and (47) into (39): 



S(m, <p; A) = 

= S harm (m,A) + 1 - lnp - 2 d ~V (l - Qm{A) 
= S(cp) + S harm (m,A) + 2 d - 1 i P g m (A) 



In p 



(1 



2 

2 d-i 



rn 



2 ip - 
- mm) 



-(1 -m) \n{2irA/d 2 



POO 

/ dye y 




' — OO 





(50) 



From this expression (that might be exact in the large 
dimension limit) it is possible to derive the phase diagram 
of the system. 



C. The equation for A and the clustering transition 

From Eq. (50) we obtain the equation for A from the 
condition f§ = 0: 

2 d tp 1-m dA 

The solution to this equation must be substituted in 
Eq. (50) to obtain the replicated entropy S(m,cp). For 
generic m, J- m {A) has the shape reported in Fig. 12 for 
m = 1. Therefore a solution exists only if 



> 



1 



2 max^J" m (A) 



<Pd(m) , 



(52) 



which defines the clustering transition density ipd{m). 
As max^ T m (^4) is a quantity of order 1 for all m, the 
scale of the clustering transition is cpa ~ d/2 d . The 
same result was obtained in (Kirkpatrick and Wolynes, 
1987a) by means of density functional theory. Moreover, 
this result compares well with the results of Torquato et 
al. (Torquato et al, 2006) who found that simple algo- 
rithms are able to construct packings up to this scale of 
density. 
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For ip > ifd(m), the equation for A admits two so- 
lutions, the physical one being the smallest, that cor- 
respond to a maximum of the free energy for to < 1, 
as usual in the replica computations. This solution has 
the right physical behavior as the cage radius becomes 
smaller on increasing the density. While the full curve 
ipd{m) can easily be computed numerically from (52), we 
will focus in the following on the special cases to = 1 and 
to = that define the equilibrium clustering transition 
(fid and tpth, respectively. 



1. The clustering (Mode-Coupling) transition 

For m = 1, the clustering transition can be identified 
with the usual dynamical (or Mode-Coupling) transition. 
This is because when the liquid splits into an exponen- 
tial number of glassy states, we expect the dynamics to 
become very slow as the system must cross barriers to 
change state. In d — > oo these barriers become very high 
and we expect a real dynamical transition characterized 
by a divergence of the relaxation time at <pd- It would be 
very interesting to check this conjecture by investigating 
the Mode-Coupling equations in the limit d —¥ oo. 

We have from Eqs. (51), (47), taking 26 the limit m -> 1 



-a^ r 

OA J- 



dye y Q 



y + A 



ln6 



y + A 



AA 



This function can be easily computed numerically: the 
result is shown in Fig. 12. It has a maximum for A = 
0.576 where T\ = 0.208; the corresponding value for the 
clustering transition is 



¥>d = ^(1) =4.8 



2 d 



d — > oo . 



(54) 



Note that this value is somehow larger (but not too much) 
than the saturation density (Prsa = d/2 d (Torquato 
et at, 2006) of the Random Sequential Addition process 
in which spheres are added sequentially at random. 



2. The J-point in d — > oo 

In the limit to — > the cage radius A that optimizes 
the free energy goes to 0. Therefore the jammed states 
are obtained in the double limit to— > 0, A 0, A ~ am, 
as will be discussed in section VI. D. Using the results of 
Appendix C, in particular Eq. (C35), one can show that 



26 Before taking the limit it is convenient to substitute in Eq. (47) 

8(y) — ¥ © ( v f A ^_ ) ; it can be shown that the integral is un- 
\V4aJ 

changed. This is due to the property f dr[qA(r) — x( r )] = 0, 
see Appendix C. 



(55) 



the equation (51) for A becomes in this limit 

d i r 00 2 

_ =7 - o(S) = _y o dyje-y-h. 

The function F^{a) has the same qualitative shape of 
^(A), see Fig. 12, and assumes its maximum for a = 
0.302 with J-o = 0.160. This gives the leading order in 
d — > oo of the clustering density for m = 0: 



d- 



oo 



(56) 



As previously discussed this density correspond to the 
first appearance of jammed states at infinite pressure and 
we conjecture that it coincides with the J-point density, 
at least in this large dimension limit. 



D. The ideal glass state 

Once Eq. (51) has been solved, one can obtain the 
entropy and complexity of glassy states from Eqs. (6) 
and (50): 



s*(m, ip) — —d\nd 



ln(27rA) + 1 + 



m 



(57) 



+ 2 d - 1 v d m g m {A) , 

Z(m,<p) = pnd-±[l + ]n(A/m)} 

- 2 d -V[i + rr?d m mT x G m {A)\ 

- ln(2 V) - \ In(Trd) + 1 , 



where we used the Stirling formula to expand the Gamma 
function appearing in Vd- In the region of the clustering 
transition we already showed that 2 d ip <~ d, and A = 
O(l); thus the first term in E always dominates and gives 
a positive complexity. 

To find the solution to S(to, ip) — we have then 
to choose f = T^jip and look to <p — > oo. In this case 

the solution for A vanishes. For small A, the function 
Q m {A) ~ V AAQo(m), see Appendix C. Substituting in 
Eq. (51) we get A oc <p~ 2 , and substituting this in (57): 

£(m,£) = ^lnd-^+0(ln£)] (58) 

Then we see that the leading behavior is (p = In d inde- 
pendently of to. The dependence on to comes from the 
subleading terms, see (Parisi and Zamponi, 2006b) for an 
analysis of these terms. Therefore, both the Kauzmann 
and glass close packing density scale as 

Vk, <Pgcp = —^d— . d -> 00 > I 59 ) 



for d — > oo. Note that in the limit to — > 0, A = am, the 
vibrational entropy goes to — oo, due to the term In A in 
(57). In summary, for to — > the cage radius goes to 0, 
the entropy goes to — oo and the pressure diverges: the 
particles get in contact and the system is jammed. 
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E. The correlation function 

The correlation function of the glass is given by the 
correlation function of the effective Baxter liquid. In the 
ring resummation for Hard Spheres, Eq. (42), we get 

I SS 



9(r) 



y(r 



dq 

Y° 



X(r)y(r) , 



(60) 



-iqr _ 



1 - Pf(q) 



where f(q) is given by (43). For the Baxter effec- 
tive liquid one has to substitute f(q) — > feffil) an d 
X(r) — > e~^ e ^ f ^ r \ j n DO th cases, for d — > oo, the analysis 
of (Parisi and Slanina, 2000) or a direct investigation of 
(60) shows that y(r) = 1 + e~2?( r ) ; where the function 
y(r) can be determined, see (Parisi and Slanina, 2000). 
Ignoring this exponentially small correction 27 we are left 
with gcir) — e~^ eff ^ r \ Therefore the correlation func- 
tion of the glass is essentially given by a step function 
plus a peak in r = D, which becomes a delta function in 
the jamming limit m — > 0, A = am. We can compute 
the weight of the delta peak x(r)Q(r), which is related 
to the average number of contacts per sphere z by: 



P / drx(r)Q(r) 



pfl d D d 
12r 



2 d <fG m (A) . (61) 



In the jamming limit Q m {A) — > Go(ce), see Eq. (C36), and 
a = a(ip) is the solution of Eq. (55), so that 



z = d 



Go(a) 



J* (S) 



(62) 



where To(a) 



- rfgo(S) 
da 



is defined in (55). In the limit 



2^0, that corresponds to random close packing, we 
have Go (a) oc Va which implies z — > 2d, i.e. the packings 
are isostatic. However, for a > Eq. (62) gives z > 2d, 
with z ~ 3.6d at the threshold density. 

This strange result might be due to different reasons: i) 
subleading corrections that we neglected might affect the 
result in some subtle way; ii) the states corresponding to 
finite a might be unstable; it is possible that only states 
with a <~ (lnd)~ 2 are stable; Hi) it is possible that other 
contributions, e.g. those of the square root singularity in 
lower dimension, are merged with the true contacts on 
the scale of the delta peak. This point deserves further 
investigation. 



F. Discussion 

We showed that for d — > oo the model can be almost 
completely solved, at least within our initial assumptions, 



i.e. that there is a clustering transition where the con- 
figuration space splits into many disconnected clusters 
without further structure 28 , and the phase diagram com- 
puted in full detail. 

To further check the consistency of the small cage ex- 
pansion, it is interesting to estimate the Lindcmann ratio 
in the glass phase, when 2 d ip ~ dlnd. The Lindemann 
ratio L for a given solid phase is the ratio between the 
typical amplitude of vibrations around the equilibrium 
positions and the mean interparticle distance. In our 
framework it can be defined as 



L = p^ d ^A . 



(63) 



Using \J~A = dy/A/D cx l/(p oc (Ind) 1 as derived in 
section VI. D, and 2 d ip = pVd(D) ~ dlnd, one has 



Vd\nd 



«1 , 



(64) 



which is consistent with the assumption that vibrations 
are very small. 

We can compare our prediction for the glass close pack- 
ing (ficcp ~ 2~ d d\nd with the best available bounds 
on the density of crystalline packings. A classical lower 
bound for the close packing density is the Minkowsky 
bound, if ~ 2~ d . It has been improved in (Ball, 1992) 
where, for the case of lattice packings, it is proved that 
if > 2d2~ d ; see also (Krivelevich et at, 2004) where a 
procedure to construct packings achieving this bound has 
been discussed. The best currently known upper bound 
is reported in (Kabatiansky and Levensthein, 1978), and 
has the asymptotic scaling <p ~ 2~°- 5990 --- d Our result for 
ifGCP lies between these bounds (it is only a factor ln<i 
better than the lower bound) so we cannot give an an- 
swer to the question whether the densest packings of hard 
spheres in large d are amorphous or crystalline. Hope- 
fully better bounds on the density of crystalline packings 
will address this question in the future. The values of 
densities of crystalline laminated lattices (Conway and 
Sloane, 1993) up to d=50 seem to suggest that there are 
lattices where 2 d ip grows exponentially in d. It is how- 
ever quite possible that this is a preasymptotic effect; 
see (Parisi, 2008) for a detailed discussion. It would be 
very interesting to find the density of laminated lattices 
in larger dimensions. Finally, it is worth to note that 
in (Torquato and Stillingcr, 2006b) it has been proposed 
that it is possible to achieve packings having density ex- 
ponentially higher than the Minkowsky lower bound, and 
actually very close to the upper bound cited above. Try- 
ing to prove (or disprove) this conjecture is a challenge 
for future research. 



That however is necessary to ensure that S(q) > for all q and 
that 5(0) = (Parisi and Slanina, 2000). 



In replica jargon this corresponds to a 1RSB solution. Unfortu- 
nately we cannot study the instability of this solution towards 
more complicated replica symmetry breaking solutions. 
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VII. FINITE DIMENSION: FIRST-ORDER IN THE 
SMALL CAGE EXPANSION 

In this section we will discuss the simplest approx- 
imation that works well in finite dimension, namely a 
first-order expansion in the cage radius \J~A. Despite its 
simplicity, we will see that it is able to reproduce many 
of the available numerical data with good accuracy. Its 
main drawback is that it does not allow to access the clus- 
tering transition, nor the full shape of g(r) in the glass 
phase. In section VIII, we will discuss different ways to 
try to improve over this approximation. 



A. First order replicated free energy 

We will focus on the first correction in the expansion 
of Eq. (40) in powers of \[~A. Then, we can neglect all 
the n > 3-body potentials that give contributions 0(A), 
see Appendix B. Moreover, at fixed m the two-body 
potential is 0(yA) and we can treat it as a perturbation. 
Note that in the limit m — > this is not true, as in fact 
the expansion parameter is A/m and not A. However we 
are interested in the limit m — > with A = am. Taking 
the limit m — > of the first order term in \[A will then 
give the first order term in ^/a. 



1. The small cage expansion 

Taking into account only the two-body potential, 
Eq. (40) reduces to (39), where 

-PF eJ f[p;<t> eJf (r)] = JV-^xWa + QW)] , (65) 

where S[ip; b(r)] is the entropy functional (32) for a non- 
replicated liquid with interaction b(r) = e^^'fJ^. It 
is easy to show, using standard liquid theory (Hansen 
and McDonald, 1986), that for a translationally invariant 
system 



1 5S[<p,b(r)] _ p 
N S\nb(r) ~~ 2 



9(r) , 



(66) 



where g(r) is the pair distribution function of the liquid. 
Using this relation to expand (65) for small Q(r), and 
calling S(<p) — N~ 1 S[(p; x( r )} t ne entropy of the non- 
replicated hard sphere liquid, the first order expression 
for the entropy is obtained from (40): 



S(m,<p;A) =7V-\S[p(x),x(x,y)] 



= S harm (m,A) + S(ip) + 



P 



J drg(r 



)Q(r) ■ 



(67) 



To compute the correlation function of the glass, we 
start from Eq. (B12); we use Eq. (67) substituting x 

Xe~ v and we get 



2 2 2 

P , \ P , ,,P 

■j9G{x,y) = -^g(x,y) + — 



, Sg(u,v) 
auav— ; -Q(u,v) , 



5hix(x,y) 



and using Eq. (D9) of Appendix D: 

6g(u,v) 



5\n X (x,y) 



g(x,y)5(x— u) 6 (y—v)+ other terms , (69) 



(where the other terms contain basically the connected 
4-point function), we get 

g G (r) = g(r)[l + Q(r)] + other terms (70) 

which is the result we already obtained in (Parisi and 
Zamponi, 2005) in a completely different way. One can 
argue that the contribution of the other terms is negli- 
gible for r <~ D; therefore Eq. (70) allows to access the 
shape of the contact peak in ga{ r )- 

if Va <c d we can assume that the function g(r) is 
essentially constant 29 on the scale V^4, g(r) ~ <?(-D)x( r ) 
for r - D - 0(\/A). Then we have 



drg(r)Q(r) ~ g(D) / dr X (r)Q(r) = g{D)V d {D)G m {A) , 

(71) 



and 



S(m,<p;A) = S harm (m,A) + S(<p) + g(D)2 d - L <pG m {A) , 

(72) 

which resembles the expression for the Baxter model in 
d — >• oo (50), except for the factor g(D) in front of the last 
term. The expansion in powers of y/~A is easily obtained, 
from Eq. (C23) in Appendix C: 

S(m, <p; A) = S harm (m, A) + S(<p) + 2 d ^g(D)^0-Q o (m) 

(73) 



This result has been obtained in (Parisi and Zamponi, 
2005, 2006a). 



2. Optimization with respect to the cage radius 

We finally have to optimize with respect to A in 
Eq. (73) to get the free energy at the minimum. Note 
that the derivative of (73) is a linear function of VA, 
therefore it always allows a solution independently of the 
values of ip and m. Therefore investigating the dynamic 
transition, i.e. the point at which the solution for A dis- 
appears, is not possible within this approximation. This 
is because the function F m (A) has a shape similar to 
the one in figure 12: if we expand it for small A and keep 
only the leading \fA term, we lose the maximum in figure 
12 and are left with a continuously increasing function. 



(68) 



29 In fact this assumption is false: it will turn out at the end of 
the computation that \J~AjD ~ 0.01 but on this scale g(r) still 
has structure close to contact. However this assumption gives 
consistently the first order in the small cage expansion. 
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Therefore the equation F m (A)=const. will have solution 
for all values of the constant, i.e. of the density. 

Keeping this in mind, we obtain for the replicated free 
energy: 



S(m, ip) = S(ip) - ^(1 - to) ln[27rA 



1 — TO 



D 



d 
2 

1 — TO 



(to — 1 
D 



In to) 



Q (m) P V d {D)g(D) Q (m) 2*<pY(<p) 



(74) 



where 30 Y(ip) = g(D). Applying Eq. (6) we get 

s*(m,<p) 



dS(m,ip) d 

= — In \2nA to 

dm 2 1 V n 



+ d(i- m) 9jM + d -^±l 

Qo{m) 2 to 
£(m, ip) = S(m, <p) — tos*(to, ip) 

= %>)-^ln[27rA*(m)] 

_ rfm (l_ m )goM + ^l nm . 
Qo(m) 2 



(75) 



dm 



The important remark is that the only input in the ex- 
pressions above is the entropy of the liquid, S((p); the 
only other quantity is Y(p>) appearing in A*(m), see 
Eq. (74), but the latter is related to S(ip) by the gen- 
eral relation 

P =^ = l + 2^ 9 (D) = -/-^ (76) 
p dip 

where p is the reduced pressure. 

Given S(<p), the density (pocp is the solution of 
Sj(<£>) = lim m _j.o S(?n, ip) = 0; using Eq. (75) and the 
asymptotic behavior of Qo(m) ~ ^/-7t/(4to) (see Ap- 
pendix C), we obtain the condition: 



S,(v?) = lim E(m, ip) = S(<p)~d\n 

m— >0 



V8 



2 d ipY{p) 



= . 



(77) 



while the Kauzmann density is the solution of Yi eq {ip) = 
S(l, if) = 0, which gives the condition 



Z eg (<p) = S(ip) - din 



2tt 



2 d Q <pY(<p) 



, 



(78) 



with Q = -Q' (m = 1) = 0.638 . . .. 

As Qa(m) ~ m -1 / 2 for to — > 0, from the second 
Eq. (74) one can show that in the limit to — > 0, ^4 cx to, 
and the first Eq. (75) shows that the vibrational entropy 



The reason for this notation is that y(r) 



Mr) 



g(r) is continuous 



m r = D (Hansen and McDonald, 1986) and Y(ip) = y(D) 
9(D). 



d 




<Pgcp 


<pj 


<Pmrj 




tpGCP 




(theory) 


(theory) 


(num.) 


(num.) 


(extr.) 


(extr.) 


2 


0.8165 


0.8745 


— 


— 


— 


— 


3 


0.6175 


0.6836 


0.640 


0.64 


— 


— 


4 


0.4319 


0.4869 


0.452 


0.46 


0.409 


0.473 


5 


0.2894 


0.3307 


— 


0.31 


— 


— 


6 


0.1883 


0.2182 




0.20 






7 


0.1194 


0.1402 










8 


0.0739 


0.0877 










oo 


2- d dlnd 


2- d dlnd 




2- d d ? 







TABLE III Values of ipx and ifiacp from the first-order small 
cage expansion of section VII (only values for d < 8 are re- 
ported for brevity, values for d > 8 are in Fig. 13) compared 
with the available numerically measured values of (fij (O'Hern 
et al, 2003; Schreck and O'Hern, 2008; Xu et al, 2005) and 
of (fiMRj (Skoge et al, 2006). The last two columns give the 
values of ifiK and (pacp extrapolated from the fits of the data 
of (Skoge et al, 2006) reported in figure 2. The large-d scaling 
of (fMRj has been conjectured in (Skoge et al, 2006). 



and the pressure both diverge. This is consistent with 
the identification of the limit m — > and the jamming 
limit, as we already discussed in section VI. D. 

It has been shown in (Parisi and Zamponi, 2005, 2006a) 
that these first order expressions already give quantita- 
tive predictions that are in very good agreement with 
numerical simulations. In the following we will review 
these results. 



B. Equation of state of the glass and complexity of the 
liquid 

Once an equation of state for the liquid is given, 
Eq. (75) give access to the phase diagram, except for 
the dynamic transition. In this section we discuss the 
results in different space dimensions. 



1. d= 1 

In dimension d = 1 we do not expect any glass transi- 
tion. Indeed, the particles are always in a "state" : apart 
from the trivial symmetries (permutations and global 
translations), they vibrate around the "configuration" 
x\ < X2 < ■ ■ ■ < xn. The packing fraction is tp = pD. 
The diameter of the cages is given by 2R = D-^, and it 
is easy to show that the distribution of the cage radii is 
exponential. Therefore in this case we expect our Gaus- 
sian approximation to be very inaccurate. 

The entropy can be computed exactly and is given by 
S((p) = 1 + ln(2i?) = 1 - \n(ip/D) + ln(l - tp). The value 
of the correlation function at contact is Y(ip) 
diverges at close packing ip = l. 

In our approach, as there is only one state, we ex- 
pect to have Z m = exp[mNS(ip)], then S(m,<p) = 



ih^ and 
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with very small quantitative differences. In the first case 
we find 31 (fx = 0.816 and <pgcp — 0.874, while in the 
second we find ifK = 0.811 and (fiacp = 0.873. 

Note that in d = 2 the existence of an ideal glass tran- 
sition has been recently the object of an intense debate 
(Brumer and Reichman, 2004; Donev et at, 2006; San- 
ten and Krauth, 2000, 2001; Tarzia, 2007). In fact, amor- 
phous packings of monodisperse two-dimensional spheres 
are particularly unstable (Brito and Wyart, 2006). For 
this reason, in the following we will focus on d > 3. 

3. d > 3 



10 



15 



d 



FIG. 13 Plot of ipK (open squares, obtained solving Eq. (78)), 
fGCP (full diamonds, Eq. (77)), and (fMRJ (full circles, nu- 
merical estimate of (Skoge et al, 2006)) as a function of the 
dimension. Both ipx and <pgcp scale as 2 d ip ~ dlnd for large 
d, while their distance scales as 2 d [ipGCP — fx] ~ d. In the 
inset the same plot for 3 < d < 6 (compare with Fig. 6 in 
(Skoge et al, 2006)). 



1 lnZ m = —mS(ip), and s*(tp) — —d m S(m,ip) — 
the vibrational entropy is equal to the total 



—N 

S(<p), i.e 

entropy and the complexity is £ = 0. Indeed, substitut- 
ing the expressions above for S(ip) and Y(tp) in (75) we 
easily find: 



2y/A*{\) = D 



1 1-ip 
0.638 ip 



a » = -ln[27rA*(l)] 



(79) 



In- 



2 • 0.638 



- Hfp/D) + ln(l - <p) 



Thus the theory reproduces the exact result apart from 
a (small) constant shift of the entropy (corresponding to 
a multiplicative factor in front of the cage radius) . This 
is probably due to the Gaussian approximation, that is 
clearly wrong as the distribution of the cage is exponen- 
tial, and to the small cage approximation. Indeed, for 
ip ru 1, where the cages are small, the leading term of the 
entropy is correctly reproduced by the expressions above. 



2. d- 



In dimension d = 2 we can use either the Henderson 
expression for Y(ip) (Henderson, 1975): 



Y H {<P) 



1 - 7^/16 



(1-^)2 

or the improved expression of Luding (Luding, 2001) 



(80) 



Y L (ip) = Y H (<p) 



•r 



27(1 - v y 



(81) 



In d = 3 we used the Carnahan-Starling expression 
(Hansen and McDonald, 1986) for the entropy S((p), 
which reproduces very well the numerical data for the 
equation of state of the hard sphere liquid. In d > 3 the 
Carnahan-Starling equation of state can be generalized 
(Song et al, 1989): 



Y(<p) = 



1 



aip 



(1 ' 

a = d-2 d - 1 {Bz/b 2 ) , 



(82) 



where Y(ip) = g(D) is the value of the radial distribution 
function at contact, and b and -B3 are the second and 
third virial coefficients, whose exact expression is known 
(Song et al., 1989). The entropy of the liquid S(ip) is 
obtained by integrating the exact expression (76). In 
d — 4 this equation is not very accurate, see figure 15, 
and an equation of state based on Pade approximants 
(Bishop and Whitlock, 2005) seems more accurate, see 
figure 1. Still the error is not so large and for simplicity 
in the analytic computations we will use the Carnahan- 
Starling approximation for all d > 3. 

Using this expression for S(<p), Eqs. (77) and (78) can 
be easy solved numerically to get the values of <Pgcp and 
fK for any given value of d. The results are reported in 
table III for d < 8, and compared with tfMRJ as reported 
in (Skoge et al., 2006). The latter quantity is the density 
of the maximally random (according to some measure of 
"order" ) collectively jammed packings of the system, see 
(Torquato et at, 2000) for the precise definition; it is 
estimated in (Skoge et al, 2006) by the jamming density 
fji'i) for finite but small 7 (see sec. IV in (Skoge et al., 
2006) for a detailed discussion). As ^(7) is expected to 
increase on decreasing 7 and ifacp = lim-y— j-o ( see 
figure 2), it follows that <pmrj, as estimated in (Skoge 
et al, 2006), should be strictly lower than (pacp, but 
close to it, consistently with the data in table III. A plot 
of ipx and ipocp f° r d up to 19 is reported in Fig. 13. 



31 Note that in Ref. (Zamponi, 2007) it was erroneously stated that 
in d = 2 the present method predicts the absence of a glass 
transition. We wish to thank F. Caltagirone for pointing out the 
error and providing the correct values of ip k and fGCP- 
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FIG. 14 Equation of state in d = 3 using the Carnahan-Starling equation of state to describe the liquid, and in the small cage 
approximation (section VII). (Left) (m,ip) phase diagram (compare with figure 6). The line m a (<p) (dot-dashed-dashed red 
line) separates the liquid and glass regions. The glass transition density is defined by m s (<PK) = 1 and the glass close packing 
density by m s (<PGCp) = 0: we find <£>k = 0.62 and ipgcp = 0.68. The blue dot-dashed line refers here to a metastable glassy 
state defined by E(m, tp) = Ej (here we chose H, = 1.2) with jamming density tpj = 0.659. (Right) Inverse reduced pressure as 
a function of density. Red full line is the CS equation of state; red dot-dashed-dashed line is the ideal glass; blue dot-dashed line 
the same metastable glass of left panel. Data from (Rintoul and Torquato, 1996) (full circles) and more recent ones from (Skoge 
et ai, 2006) (dashed black line, same data as in figure 1, left, for 7 = 64 ■ 10 -6 ) are reported for comparison. 



Note that it has been suggested in (Skogc et at, 2006) 
that 2 d (fMRj ~ c\ + c%d] and indeed, given that the 
compression rates used in numerical simulation are not 
very small, (fMRJ should be quite close to ifth, that scales 
as d/2 d , see Eq. (56). Recalling that 2 d ip K ,2 d ip GC p ~ 
dlnd, we expect that at some point (Pmrj will become 
smaller than tpx, even if this is not observed for d < 6. 

The very nice data for d = 4 reported in (Skoge et at, 
2006), and reproduced in figure 1, allow for a more pre- 
cise comparison of the numerical and theoretical results: 
the value of ^(7) has been measured for 5 different val- 
ues of 7 = 10~ 3 , 10~ 4 , 10- 5 , 10~ 6 , 10~ 7 , see figure 2. A 
standard procedure to extrapolate to 7 — > is to fit the 
data to a Vogel-Fulcher law: 

o D 

= 70l0 "GCP-Vi & y>.( 7 ) = tp GC p + - —— . 

Such extrapolations are often ambiguous; however the fit 
is good (see figure 2) and gives <pgcp = 0.473. A simi- 
lar extrapolation of tp g (7) (defined roughly as the point 
where the curves in figure 1 leave the liquid equation of 
state) to 7 = gives <px = 0.409. The final results dif- 
fer by ~ 10% from the theoretical values, see table III. 
This is a very good result given the ambiguities that are 
present both in the numerical data (numerical error and 
extrapolation) and in the theory (the choice of a par- 
ticular expression for the equation of state of the liquid 
that is not exact, and the small cage expansion). Note 
that a similar extrapolation is not possible in d < 4 due 
to crystallization, and for d > 4 due to lack of numeri- 
cal data. Hopefully new data for d > 4 will allow for a 
similar comparison also in this case. 



In figures 14 and 15 we report the inverse reduced pres- 
sure as function of the packing fraction for d = 3 and 
d = 4. We plot the numerical data of (Skoge et ai, 2006) 
already reproduced in figure 1, and we choose two among 
the smallest compression rates available. Still we see that 
the numerical data are far from the ideal glass pressure, 
and are better described by metastable glassy states cor- 
responding to a complexity £ ~ 1. Indeed, it is well 
known from numerical simulations of structural glass for- 
mers that the system falls out of equilibrium when E ~ 1. 
The overall agreement of the theoretical predictions with 
the data of (Skoge et ai, 2006) is very good. 

Finally, in figure 16 we report the equilibrium complex- 
ity H eq (ip) = E(l,<£>) as function of (p. Numerical data 
from (Angelani and Foffi, 2007) and (Speedy, 1998) are 
available, and the agreement is again very good. 



C. Scaling close to jamming 

In this section we will derive asymptotic relations for 
the behavior of the glassy states at ip — > ipj . This is inter- 
esting because close to ipj the correlation function shows 
interesting features that have been studied in detail nu- 
merically, see the discussion in section II. D. 



1. Scaling of the pressure 

First we compute the scaling of the pressure close to 
ipj. It is convenient to recall the procedure already dis- 
cussed in section III. A. 3. To each group of glassy states 
of jamming density ipj we can associate a complexity Ej 
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FIG. 15 Inverse reduced pressure as a function of density in 
d = 4 using the Carnahan-Starling equation of state (red full 
line) to describe the liquid, and in the small cage approxi- 
mation (section VII). Data from (Skoge et al, 2006) (black 
dashed line, same data as in figure 1, right, for 7 = 10~ 6 ) 
and from (van Meel et al, 2009) are reported for comparison. 
The reported metastable glass (dot-dashed blue line) corre- 
sponds to Ej = 1.44 and has a jamming density ipj = 0.467. 
In this case the Carnahan-Starling equation is less accurate; 
nevertheless, the overall quantitative agreement is still good. 
The ideal glass branch is reported as a dot-dashed-dashed red 
line. 



given by Eq. (77). Then one has to solve Eq. (9) to 
obtain m(tp, Sj) or equivalently m(tp,tpj). The function 
m(tp, ifij) represents a group of glassy states in the (to, tp) 
plane of figure 6, left panel. Once substituted in equa- 
tion (75), it gives the entropy of the states 32 labeled by 
tpj as a function of tp. 

We are interested in the equation of state of these 
states close to m = 0, or tp = tpj. Then we need to 
compute m(tp,tpj) at first order in tpj ~ tp. To this aim 
we have to linearize S(m, tp) — Tij{tpj) at first order in m 
and tfj — tp; after a tedious computation (see Appendix E) 
we get 



T,(m,tp) - Hjifj) 



(<Pj -<p)-dm , 



(83) 



that gives, expressing S (tp) in terms of Y(tp) by means 
of Eq. (76) 



m(tp,tpj) 



,» 1 „ s Y'ftpA l-d 
2 d - 1 Y(tp 3 ) - d — ^ 1 



Y(tpj) tpj 



(<Pj - <P) 
(84) 
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0.46 0.48 0.5 0.52 0.54 0.56 0.58 0.6 0.62 

FIG. 16 Equilibrium complexity E e9 (<p) as function of tp in 
d — 3 using the Carnahan-Starling equation of state in the 
small cage expansion (section VII). The prediction is com- 
pared with data from Angelani and Fofli (Angelani and Foffi, 
2007) and Speedy (Speedy, 1998). 



The latter result must be substituted in s* (m, tp) to get 
the behavior close to ipj. It is easy to see from Eqs. (75) 
that for m — > 0, s*(m, tp) ~ d(l — m) In to; therefore 



s(tp, ipj) - dln(tpj -tp)- dn(tpj)(tpj - tp) ln(tpj -tp) + - 



close to tpj, and the (reduced) pressure diverges as 

ds(tp,<pj) 



dtpj 



dtp 



(85) 



(86) 



dtpjKv^Mvj -v) + 



It is important to remark that the corrections to the lead- 
ing order are quite large (logarithmically divergent), as 
observed in numerical simulations (Skoge et a/., 2006). 



2. Correlation function 

The most interesting results on the scaling close to tpj 
are obtained by investigating the pair correlation func- 
tion of the spheres in the glass state. At first order we 
obtained Eq. (70); the first term describes the delta peak, 
while the other terms do not contribute close to r = D. 
It would be very interesting to evaluate them, but this 
seems a very complicated task (see Appendix D for a 
complete list of these terms). Here we will focus on the 
first term. We have then to study the asymptotic be- 
havior of Q(r) in the jamming limit. This is done in 
Appendix C.3, and the result is 



g G (r)=g(r)[l + Q(r)]=g(r) 



1 



1 



-A„ 



r-D 



32 We recall that the case Sj = correspond to the ideal glass. 



where g(r) is the (extrapolated) correlation function of 
the liquid at the same density, and to = m(tp,tpj), A — 
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and therefore 



l + 2 d - 1 V9G {D) 



■jd-l 



dipj 



o.i 



10 



p(<p,(pj) . 

(89) 

This is not a surprise since the relation (76) is vio- 
lated in most of the approximate theories for the liq- 
uid phase (Hansen and McDonald, 1986). Note however 
that the difference is not so big, because the first term 
in (/.((fj), given by Eq. (84), dominates as Y(ip) is quite 
large (~ 20) in the density range of interest. 

Using again m(cp, ipj) = fj,(ipj)(tpj — ip) and the expres- 
sion of A*(m) we get from Eq. (87) close to <py. 



FIG. 17 Scaling limit of the delta peak. Numerical data, 
reproduced from (Donev et al., 2005a), refer to a packing 
(fj ~ 0.64 and <Pj — <p ~ 10~ 12 . ga(r)/gG(D) is plotted against 
A = ^-jy-p. The full line is the theoretical prediction (90). 
Compare with figure 9 for a similar result obtained in the 
HNC approximation. 



A* (m(ip, ifij)); the function A m (t) is defined in Eq. (C38) 
and plotted in figure 25. 

The function A m (t) gives the leading contribution to 
9g{ t ) close to contact, that comes from neighboring par- 
ticles that are in contact with the reference particle ex- 
actly at ifj. This contribution is nonvanishing only for 
r — D < y/A ~ y/ifj — ip, and becomes a delta peak at (fj . 
Other corrections are encoded in the remaining terms in 
(70) or are missed by the first order small cage expan- 
sion. A first consequence of this fact is that the integral 
°f 9G( r ) ~ 1 does not vanish at <p = ipj as it should, 
being proportional to the compressibility. A second con- 
sequence of our approximations is that Eq. (76) is not 
satisfied: this is because from Eq. (87), for <p — > ipj 



g G (D)=Y(<p)[l + Q(D)] 



npj) 



9a{r) 
9g(D) 



(r - D^QtpjYfa) 



y/n r — D 



V) 



D 



{l + 2 d - x V g G {D))\ (90) 



Aq 
A 
Aq 



having defined the variable A = ^(1 + 2 d ~ 1 >pg G {D)) ~ 
^-jp-p (the last equality holds only if (76) holds). Note 
that the same scaling of the delta peak was already pre- 
dicted in the HNC approximation, see figure 9; thus it 
is a very robust feature that emerges from the replica 
method independently of the approximations. 

This scaling form of the delta peak of g G (r) close to 
jamming has been observed numerically in (Donev et al, 
2005a). In Fig. 17 we report numerical data from (Donev 
et at, 2005a) on the delta-peak contribution for a packing 
in d = 3 with ipj ~ 0.64 and <pj — <p~ 10~ 12 . We see that 
our result, given by Eq. (90) for ipj = 0.64, is in excellent 
agreement with the numerical data. 

In summary, from Eq. (87) and Appendix C.3, we get 
the following scaling for the delta-peak of g G (r): 



9G{r) 



k exp[-(r- Df/(ipj - ip)} 



r — D ipj — ip , A~ 1 , 

<Pj-ip<r-D~ y/pj - ip , A - l/y/ipj - ip 



(91) 



r - D » y/ipj - <p , A > 1/y/pj - ip , 



where the last line comes from the cutoff on the function 
A m (t) at finite m, see Appendix C.3. 

In numerical simulations a (r — D)~ a divergence of 
g G (r) close to D is observed, with a ~ 0.5. This leads 
to a difference with respect to (91), with the (r — D)~ 2 
regime extending only up to r — D ~ (ipj — p) 2 ^ 3 . This 



feature is not present in our calculation; this discrepancy 
deserves further investigation. 



3G 




FIG. 18 Number of contacts from Eq. (95). (Left) Data from Ref. (Donev et al, 2005a) in d = 3 and tpj ~ 0.64. We have 
z ~ 6 + 20yVj — tp for both curves. (Right) Data from Ref. (Brito and Wyart, 2006) in d — 2 and tpj ~ 0.83. The best fit gives 
z — 4 + — tp while the theory gives z = 4 + 14.5^/tpj — tp. 



3. Number of contacts 

The interpretation of Eq. (91) is that for tp — > tpj there 
is a shell of width tpj — p around a given particle where the 
probability of finding other particles is very high. These 
particles become neighbors of the particle in the origin 
at p — ipj, and their number is finite as the integral of 
9g{t) on this shell is finite. 

There is however a second shell tpj — tp < r — D < 
^ifj — ip where the probability of finding particles is very 
high. The integral of go( r ) over this shell is of order 1 so 
also in this shell there is a finite number of particles that 
will become neighbors of the particle in the origin. 

The integral of gaif) on a shell D < r < D + 
0(yj(fj — tp) gives then the number of contacts for ip — > 
iff 
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dr [1 + Q(r)] 



(92) 



Recall that in the expression above, as everywhere in the 
paper, g{D) — Y(p) is the contact value of the liquid 
correlation, that is finite; this should not be confused 
with ga{D) which diverges at jamming. For ip$ — if. finite 
but small this number can be interpreted as the number 
of particles that collide with the particle in the origin 
during a finite but long time r (Brito and Wyart, 2006). 

Remarkably, the integral of the second term in Eq. (92) 
can be computed exactly using Eq. (C39); since the func- 
tion A m (t) contains an exponential cutoff at r — D ~ \/~A 
(see figure 25), the upper integration limit can be ex- 
tended up to oo and we obtain 

POO 

n dP D d - 1 g{D) drQ(r) = 2d(l-m) (93) 



JD 

i.e. z = 2d for tp = ipj. Close to ifj there is a correction 
Sz tx m ~ ipj — p coming from the equation above, plus a 



second correction Sz cx y/tpj — tp coming from the integral 
of 1 in Eq. (92). The second correction dominates, then 
we have 



z = 2d + 0{ytp~^p) 



(94) 



as found in (Brito and Wyart, 2006). We can try to give 
a quantitative estimate of the coefficient by observing 
that the function A m (t) starts to drop exponentially at 
t = y/m and drops to for t = T^pm with r very close 
to 2. Clearly the upper limit of integration r plays the 
role of a fitting parameter, however the value r = 2 is 
very reasonable from Fig. 25, given also the uncertainty 
involved in the numerical estimate. Thus we assume that 
the integral is done up to r = D + 2y/A and we have 



z = 2d + n dP D d - 1 g(D)2VA 
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(95) 
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The result is compared with numerical data in Fig. 18. 
The value of [i(<fj) is given in Eq. (84). Despite the good 
agreement, note that the square root growth of Sz was 
attributed in (Brito and Wyart, 2006) to the square root 
singularity of go (f) , which is a different mechanism from 
the one we discussed above. 



4. Force distribution 

In (Donev et al., 2005a) the force distribution of 
dense amorphous packings generated by using the 
Lubachevsky-Stillinger algorithm in d — 3 was investi- 
gated. The interparticle forces are defined in this way: 
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one takes a packing of density (pj, slightly reduces the 
density and measures the average momentum exchanged 
by two neighboring particles over a large time t. Note 
that on very large times and if the volume is large enough 
the packing will be unstable and relax toward a more 
compact structure (either at (fGCP or Vfcc), se e 
Fig. 5 in (Donev et at, 2005a). However before this hap- 
pens there is enough time to measure the forces with 
sufficient accuracy. 

The interparticle forces are then normalized such that 
(/) = 1 and the distribution P(f) is investigated. The 
theory of (Donev et at, 2005a) relates this quantity to 
the shape of the delta peak close to contact, Fig. 17, by 
the relation 



9gW 
9g(D) 



dffP(f)> 



-A/ 



(96) 



where A = (r/D — l)p as in the previous sections. The 
previous equation expresses the fact that for small inter- 
particle separation the gap h oc 1//, as discussed also in 
(Brito and Wyart, 2006). Eqs. (90) and (C40) give: 



and finally 
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(98) 



In Fig. 19 we can see that this form reproduces well the 
numerical data for large forces. The discrepancy at small 
forces can be explained by i) the fact that the small / 
behavior of P(f) is related to the large A behavior of 
<7g(A), and as discussed above in the large A region there 
are corrections to <7g(A) that we are missing, and it) a 
possible finite size effect (indeed it seems that the data for 
larger samples are in better agreement with the theory). 

We wish to stress again that here we focused only on 
the results of (Donev et ai, 2005a) that refer to packings 
produced using slow compressions; these should be re- 
lated to the infinite pressure glassy states that are the 
object of this paper. Interestingly, similar results for 
P(f) have been obtained in (Snoeijer et at, 2004a, b) by 
using an ensemble approach for the force network based 
on Gaussian random matrices. In (O'Hern et ai, 2002) 
a Gaussian tail in P(f) was found for packings close to 
the J-point; moreover, it has been shown that finite size 
fluctuations of the jamming density can introduce self- 
averaging problems when averaging over many configu- 
rations that in turn produce exponential tails in P(f). 
Therefore some care should be taken when comparing our 
prediction of a Gaussian shape of P(f) with numerical 
data. Moreover, the force distribution has been studied 
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FIG. 19 Probability distribution of the forces, from (Donev 
et at, 2005a). The full red line correspond to the Gaussian 
form (98). The black line is a fit performed by the authors 
of (Donev et al, 2005a). 



in great detail for packings of frictional spheres produced 
using different protocols and seems to depend strongly 
on the protocol used. As discussed in the introduction, 
covering the relevant literature here is impossible. In par- 
ticular friction strongly affects the force network and is 
expected to have a dramatic effect on P(f). 



VIII. BEYOND FIRST ORDER IN THE SMALL CAGE 
EXPANSION 

In order to improve over the first order expansion, one 
could try to perform a second order (and third, etc.) 
expansion systematically. However, the convergence of 
this procedure is not clear, and moreover, the calculation 
of higher order corrections becomes increasingly difficult 
and is already very difficult at second order. Following 
the tradiction of theoretical physics, we then try to use re- 
summation techniques for our expansion. Unfortunately, 
the results presented here are not satisfactory, but still 
we believe that they are interesting since they illustrate 
a possible direction one can follow to improve the theory. 
We believe that a more accurate treatment of the resum- 
mation we will discuss could indeed lead to much better 
results. The reader interested only in "stable" results is 
then encouraged to skip this section. 

The simple resummation we will discuss in this section 
has the advantage of being fully analytically solvable, but 
unfortunately gives inconsistent results for the thermody- 
namics of the ideal glass. The main idea is to fully exploit 
the resummation discussed in section V, that has been 
particularly useful to solve the model in the limit d — > oo. 
Note that the small cage expansion of section VII does 
not give access to the pair correlation function of glass, 
which is one of the most interesting quantities to com- 
pute, see the discussion in section II. D. The advantage 
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of the resummation we discuss here is that it allows, by 
finding a good approximation to compute the properties 
of the effective liquid, a full computation of the correla- 
tion function of the glass. 



the same density ip and of interaction strength r = 
\/{AG m {A)). 



1. Free energy and the equation for A 



A. Percus-Yevick approximation for the Baxter model 

To resum a class of terms in the small cage expansion 
we start from Eq. (39), i.e. we neglect all the n > 3-body 
potentials but try to treat the two-body one as exactly 
as possible. Then we have to compute the entropy of a 
liquid with an hard core repulsion and a small tail given 
by (f> e ff(r) = — ln[l + Q(r)]. The tail is attractive for 
m < 1 and repulsive for to > 1, see Fig. If. 

The advantage of this formulation is that the correla- 
tion function of the glass is simply 
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2 5S[p, X (r)e 
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-"W(l + Q(r))] 



dv(r) 



v=0 



9eff(r) , 
(99) 

where g e ff(r) is the pair correlation of the liquid with 
the effective potential (p e ff- Therefore gaif), the correla- 
tion of the glass, turns out to be equal to the correlation 
function of the effective liquid and then automatically 
has many of the reasonable properties of a correlation 
function for hard particles (e.g. it is positive, it vanishes 
inside the hard core, the structure factor is positive, etc.). 
We can use different approximation schemes (PY, HNC, 
etc.) to compute the free energy and correlation function 
of the effective liquid. 

To further simplify the problem, we observe that in the 
limit m — > the cage radius vanishes and the strength of 
the attractive part of the potential diverges. Then we can 
approximate Q(r) with a delta function as in (45) and the 
effective liquid becomes a Baxter model. Note that the 
approximation (45) strictly holds only for to — > 0, where 
A — > 0. However it might be reasonable for any to < 1, 
as in this case the interaction is attractive and its range 
is anyway much smaller than the hard core diameter. We 
run into problems, however, for to > 1, as in this case the 
interaction is repulsive and cannot be described by (45). 
Therefore in the following discussion we always assume 
that to < f . 33 

In the notation of Baxter (Baxter, 1968) we have: 



= x(r)[l + Q(r)]->x(r) 



(100) 

Comparing this with (C4), we see that in d = 3, the 
effective liquid is finally reduced to a Baxter model at 



Remarkably, the Percus-Yevick approximation for the 
Baxter model has been solved exactly in d — 3 (Baxter, 
1968): an analytical expression for the free energy has 
been given in Eqs. (2.5) and (2.7) of Ref. (Tejero and 
Baus, f993), while the correlation function can be com- 
puted using the method of (Wcrthcim, 1963). 

The function t(to, A) can be computed in d = 3 by 
a numerical integration of Eq. (C4), (CI7). Then we 
get an analytical (although complicated) expression for 
>S(m, p; A): 



S(m,p;A) = S harm (m,A) + l-ln — -ip 

7T 



' 4G m (A) 

(101) 

where ip[p,r] is given in Eq. (2.7) of (Tejero and Baus, 
f993). The equation for A reads, using the definition of 
Appendix C.l: 



1 = 
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' 4G m (A) 



F m (A) , 



(102) 



and it has to be solved numerically for generic to, p. 

However some simplifications are possible. First of all, 
for to — > 1 we have G m {A) — > and r — > oo. It is possible 
to show that 



Spy{p) = 1 — In 
dip[(p,r] 

lim a -1 = 



6p 



lim ip[p,r] 



(103) 



-pYpy(p) ■ 



Therefore the equation for A becomes 
l = ^pY PY (p)F 1 (A) . 



(104) 



In the jamming limit to — > 0, A = am, G m (A) — > 
Go (a) that can be computed easily, see Eq. (C34). Using 
Eq. (C35) the equation for a is 



8 dip 



' 4G (a) 



F {a) 



(105) 



The expression for the complexity can be simplified 
similarly. 



2. Results 



33 Anyway, for m > 1 the strength of the interaction is always finite 
(and much smaller than for m < 1). In this case it seems very 
reasonable to use the first order approximation derived in the 
previous sections. 



Unfortunately, the results of the Baxter resummation 
are very poor in d = 3. The main problem is that the 
static value of to, m s (p) is not a decreasing function of 
p. This is inconsistent since it would imply that the ideal 
glass state at has smaller density than the liquid state at 
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FIG. 20 (Left) Phase diagram of the Baxter resummation (section VIII). The dashed (blue) line is the boundary of the 
region where a solution for A is found. The full (black) and dot-dashed (red) lines are the boundary of the region where the 
complexity is positive. (Right) Number of contacts z as a function of density along the line rrid(<p)- The black dot corresponds 
to md(ip) = 0. 2 increases almost linearly from (liquid state) to z ~ 6, the value expected for an amorphous jammed phase. 
Note however that the upper part of the line correspond to negative complexity. 



<p = <pK- On the contrary, the behavior of the line md(ip) 
is reasonable even if the value of the transition at m = 1 , 
(fd = 0.418, is exceedingly smaller than the conjectured 
one, (fid ~ 0.58. For these reasons the phase diagram in 
this approximation seems unreliable, see Fig. 20. 

More interesting is the result for the pair correlation 
function. Indeed, the pair correlation function of the 
Baxter liquid shows many characteristic features which 
are observed in jammed packings of hard spheres, as it 
has been discussed in (Miller and Frenkel, 2004b). The 
most interesting among them are the peak at r/D = V3 
and the jump in r/D = 2 which have been observed 
in (Donev et at, 2005a) for random jammed packings 
of monodisperse spheres. Unfortunately, the values of 
t which comes from the analysis above seems to be too 
high for these features to be present in the ga( r ) with the 
correct order of magnitude so that a quantitative compar- 
ison with numerical simulations seems to be impossible 
at this stage. 

Note also that the g(r) of the Baxter liquid is char- 
acterized by a delta peak at contact (due to the adhe- 
sive potential). From the amplitude of the latter peak, 
which is analytically known in the Percus-Yevick solution 
(Baxter, 1968), we can compute a number of neighbors 
as a function of density. The result is similar to what 
expected: the number of neighbors increases from to 
z ~ 6 = 2c?, i.e. an isostatic packing, see Fig. 20. How- 
ever the upper part of the curve corresponds to negative 
complexity so that the corresponding packings do not 
exist. 



B. Discussion 

This particular very simple resummation scheme seems 
not able to describe correctly the glassy states of hard 



spheres. However, the idea of resumming all the terms 
corresponding to the two-body interaction seems promis- 
ing, as the resulting gai 1 ") shows some features that are 
observed in real packings and the behavior of the delta 
peak (number of contacts) is correct. 

It is possible that three (and more) body interactions 
play a major role, however it would be worth to try to 
keep only the simplest two-body interactions and use 
more refined equations for the effective liquid. For in- 
stance, the delta-function approximation for Q(r) might 
be too rough as \f~A ~ 0.01 and on this scale the g(r) of 
hard spheres has strong variations at least close to con- 
tact. One could then try to solve numerically the Percus- 
Yevick approximation for the exact potential <f> e ff- Al- 
ternatively, a more precise equation of state (such as the 
one proposed in (Miller and Frenkel, 2004a)) could be 
used for the Baxter model. 

One might be worried because of the fact that in the 
limit m — ¥ we do not recover the small cage expansion 
results (e.g. the complexity was positive for ip < 0.68 at 
first order in \f~A while here we always found a negative 
E). The reason is that for m — > 0, the small cage param- 
eter is a — A/m and not A. Therefore, even if A — > 
in this limit, the small cage expansion is not recovered 
because a remains relatively big. 



IX. BINARY MIXTURES 

In this section we report the recent extension of the 
theory to binary mixtures that was obtained in (Biazzo 
et at, 2009). The extension to binary system is impor- 
tant since in these systems crystallization can be easily 
avoided; this allowed to study their glass transition in 
great detail, see e.g. (Donev et ai, 2006; Foffi et at, 2003, 
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FIG. 21 Inverse reduced pressure, p/(/3P) as a function of 
the packing fraction <p for a mixture with r = 1.4 and x = 1. 
Numerical data are obtained using two different protocols. 
In the first, compression is started at low density. In the sec- 
ond, compression is started from an equilibrated configuration 
at <p = 0.58. The equation of state of different metastable 
glasses, corresponding to different Ej, are reported as dashed 
lines. The dot-dashed line is the pressure of the ideal glass, 
corresponding to Ej = 0. A numerical estimate (Berthier 
and Witten, 2009b) of the Kauzmann pressure, pk = 34.4, is 
reported as a dotted horizontal line. 



2004; Gotze and Voigtmann, 2003; Santen and Krauth, 
2000) and in particular the recent works (Berthier and 
Witten, 2009a, b; Hermes and Dijkstra, 2009). Jamming 
of binary mixtures has also been extensively investigated, 
see e.g. (Dodds, 1975, 1980; O'Hern et al, 2002; Ouch- 
iyama and Tanaka, 1981). Their investigation allows to 
test the prediction of the theory concerning the varia- 
tion with composition of packing density, structure, etc. 
The details of the computation of the function S(m, tp) 
for a general multicomponent mixture, following (Coluzzi 
et al., 1999), can be found in (Biazzo et al., 2009) in 
the framework of the first-order small cage approxima- 
tion. Here we only report the results for a binary mix- 
ture of two types of three-dimensional spheres \i — A,B 
in a volume V , with different diameter and den- 
sity p M = N^/V. We define r = D A /D B > 1 the 
diameter ratio and x = Na/Nb the concentration ra- 
tio; <p = paVz(Da) + PbVs(Db) the packing fraction; 
V = PbVs(Db)/p = 1/(1 + xr 3 ) the volume fraction of 
the small (B) component. 

As in the monodisperse case, once an equation of state 
for the liquid has been chosen, one can obtain for a given 
Y>j the equation of state of the corresponding metastable 
glass, as well as its jamming density <pj. The equation 
of state used in (Biazzo et al, 2009) is a generaliza- 
tion of the Carnahan-Starling equation. The theoretical 
results were compared with numerical results obtained 
from the Lubachevsky-Stillinger algorithm discussed in 
section II. A. In figure 21 the evolution of the inverse 
reduced pressure during compression is reported for a 



FIG. 22 Packing fraction ipj as a function of T] = 1/(1 + xr 3 ) 
at fixed r. Full symbols are numerical data from (Biazzo et al., 
2009). Open symbols are experimental results from (Yerazu- 
nis et al., 1965). Lines are predictions from theory, obtained 
fixing Ej = 1.7. Note that the large r-small r\ region can- 
not be explored, since for such very asymmetric mixtures the 
large spheres form a rigid structure while small spheres are 
able to move through the pores and are not jammed (Dodds, 
1980; Ouchiyama and Tanaka, 1981). 



mixture with x = 1 and r = 1.4, that has been recently 
studied in great detail (Berthier and Witten, 2009a, b; 
O'Hern et al., 2002). Overall, one observes the same be- 
havior already discussed in the monodisperse case. The 
curve [Numerical 1) corresponds to a relatively fast com- 
pression (7 = 10~ 2 ). The curve (Numerical 2) has been 
obtained starting the compression (at the same compres- 
sion rate) from a carefully equilibrated liquid configu- 
ration of the same mixture at ip = 0.58 (see (Berthier 
and Witten, 2009b) for details on how this configuration 
was produced and equilibration was checked) . In the lat- 
ter case, since the relaxation time of the liquid at that 
density is already very long compared to the compres- 
sion rate, the system falls immediately out of equilibrium 
and the pressure increases fast until jamming occurs at 
a higher density compared to the previous case. The 
numerical equation of state is compared with that of a 
glassy states corresponding to £j = 0.5, 1.2, 1.5. Starting 
the compression from a high-density liquid equilibrium 
configuration produces a glassy state with lower Ej . This 
is a nice confirmation of a prediction of the theory, that 
different glassy states jam at different density. Finally, 
we report, for the same system, the numerically extrap- 
olated value of the reduced pressure p = f3P/ p at the 
ideal glass transition, px = 34.4, obtained in (Berthier 
and Witten, 2009b) . Again, this corresponds well (within 
10%) to the computed value px — 31.8 from the theory, 
see figure 21. 

In figure 22, the jamming density ipj is reported for dif- 
ferent mixtures, putting together numerical results (Bi- 
azzo et al, 2009), experimental data (Yerazunis et al, 
1965), and the theoretical results. The latter have been 
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FIG. 23 Partial average coordination numbers (small-small, small-large, large-small, large-large) as a function of volume fraction 
of the small particles 77 = 1 / (1 + xr [i ) for different values of r. Full symbols are numerical data from (Biazzo et ai, 2009). Open 
symbols are experimental data from (Pinson et ai, 1998). Note that in the lower right panel a different scale is used for zi s . 



obtained by fixing Sj = 1.7 for all mixtures, which gives 
the best agreement. Note that a single "fitting" pa- 
rameter Yjj, that is strongly constrained, allows to de- 
scribe different sets of independent numerical and exper- 
imental data. The prediction of the theory are quali- 
tatively similar to previous ones (Dodds, 1980; Ouch- 
iyama and Tanaka, 1981), but the quantitative agree- 
ment is much better. Interestingly, a similar qualitative 
behavior for (fMCT has been predicted by Mode-Coupling 
theory (Foffi et al, 2003; Gotze and Voigtmann, 2003); 
although there is no a priori reason why the variation 
with mixture composition of (p mct and tpj should be re- 
lated, it is reasonable to expect that they show similar 
trends (Foffi et al, 2003). 

The average coordination numbers at ipj are denoted 
z nv{fj)i and it was checked in (Biazzo et al., 2009) that 
their variations with ifj are negligible. They are reported 
in figure 23 for different mixtures. The numerical values 
have been obtained in (Biazzo et ai, 2009) by remov- 
ing the rattlers from the packing. Experimental data 
from (Pinson et al., 1998) are also reported in the right 
panel of figure 23. As discussed above, the total coor- 
dination is close to the isostatic value z — 6, which is 
the value predicted by the theory also for binary mix- 
tures (Biazzo et ai, 2009). As it can be seen from fig- 
ure 23, the computed values agree very well with the 
outcome of the numerical simulation, at least for r not 
too large, while some discrepancies are observed in the 
contacts of the large particles for large r. 



X. CONCLUSIONS 

This paper is based on the main assumption that amor- 
phous jammed packings of hard spheres can be identified 
with the infinite pressure limit of glassy states. In ad- 
dition we assumed that the mean field scenario for the 
glass transition holds also in finite dimension, at least on 
the time and length scales that are currently investigated 
in numerical simulations (and sometimes also in experi- 
ment on colloids and granular systems when the number 
of particles is not so large) (Bouchaud and Biroli, 2004) . 

The mean field picture leads to a non-trivial structure 
of the phase diagram, whose main consequence is the ex- 
istence of amorphous jammed packings in a range of den- 
sities [<fith, <Pgcp}- The random close packing density can 
be any density within this interval and its precise value 
depends on the details of the protocol used to construct 
the packings (Krzakala and Kurchan, 2007). 

Based on these assumptions we used the replica 
method to compute the properties of these glassy states. 
Before concluding we will briefly summarize our results 
and discuss some possible future developments, including 
how the picture is modified in finite dimensional systems. 



A. Summary of our results 

We used different approximation schemes for the repli- 
cated liquid. The HNC approximation seems to work 
well in the moderately dense phase close to ipj, and ipx 
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at m — 1. On the contrary, the small cage approxima- 
tion works better in the regime where cages are small, 
namely for m ~ close to the ideal glass line and in 
large dimension. In dimension d = 3 the mapping of the 
replicated liquid onto the Baxter adhesive hard sphere 
model seems a promising way to obtain a satisfactory 
description in all the phase diagram but for the moment 
gives poor quantitative results. 

We list here the main results we discussed in this paper; 
many of them have been compared with numerical results 
in the figures. 

• We presented a consistent description of the glass 
transition for hard spheres in d — > oo: in particular 
we give predictions for the clustering -Eqs. (54) and 
(56)- and Kauzmann -Eq. (59)- densities. We are 
able to compute the correlation function, Eq. (60), 
and the number of contacts, Eq. (62), which we find 
equal to 2d at least close to <pgcp- 

• We computed the Kauzmann and glass close pack- 
ing densities in any finite dimension in the small 
cage expansion, see table III and figure 13. 

• In d = 3 we obtain an expression, Eq. (97), for the 
scaling of the contact peak of g(r) close to jamming. 
This expression describe very well the numerical re- 
sults, see figure 17. From this expression it follows 
that the distribution of contact forces is Gaussian, 
see Eq. (98) and figure 19, and that the number of 
neighbors is z = 2d at jamming (i.e. packings are 
isostatic), see Eq. (95) and figure 18. The form of 
the scaling function seems particularly robust as it 
is found both in the small cage expansion and in 
the HNC approximation. 

• We are able to reproduce the equation of state of 
the glass for slow compression rates, see figures 14 
and 15, and the equilibrium complexity of the liq- 
uid, see figure 16. 

• In the HNC approximation we have indications for 
the development in g(r) of a jump in r = 2D and 
of long range correlations; these features cannot be 
studied within the small cage approximation but 
should be present also in the Baxter resummation. 

• For binary mixtures, we showed that the theory 
correctly predicts the variation with mixture com- 
position of the jamming density and the partial co- 
ordinations. 

Finally, very recently strong evidence has been reported 
for a divergence of the equilibrium relaxation time of an 
hard sphere liquid at a density ipo at which the pressure 
of the liquid stays finite (Berthier and Witten, 2009a, b; 
Brambilla et al., 2009). This seems consistent with the 
existence of a Kauzmann transition at ipo = tpx of the 
kind discussed in this paper, even if a computation of the 
complexity for this system is still missing. Still these re- 
sults seem to exclude the possibility that the equilibrium 



relaxation time diverges only at random close packing as 
proposed by some authors. 

B. Discussion and perspectives 

There are still many important points that deserve in- 
vestigation. A tentative list is the following: 

• Our theory is based on an "equilibrium" computa- 
tion (where we also include metastable long lived 
states). Hence, a key role is played by entropy. 
We wish to stress that in the context of metastable 
glassy states, there are two very important concepts 
related to entropy: the first one is the complex- 
ity, that counts the number of such states; the sec- 
ond is the surface tension between different amor- 
phous states, which must also be of entropic na- 
ture since we are dealing with hard spheres. The 
existence of a surface tension (a free energy cost) 
to match different glassy states is necessary for 
those states to be well defined; see Appendix A 
and (Bouchaud and Biroli, 2004; Cavagna, 2009; 
Xia and Wolynes, 2001a) for a more complete dis- 
cussion. While the complexity can be computed, 
at least approximately, within the theory presented 
here, a first principle method to compute the sur- 
face tension is still missing, despite preliminary at- 
tempts in simplified models (Franz, 2005). How- 
ever, the study of the latter has very recently be- 
come the subject of intense numerical effort (Biroli 
et al., 2008; Cammarota et at, 2009), which will 
hopefully lead to a more complete understanding 
of the properties of glassy states. 

• Clearly, a major open problem is whether a ther- 
modynamic glass transition really exists in finite di- 
mensional systems. Our mean-field theory has not 
much to say on this problem. It might seem that in 
our theory the existence of a Kauzmann transition 
is a major assumption. This is indeed not true: 
going carefully over the discussion, it will become 
clear that what really matters for us is the existence 
of long-lived metastable glassy states. The glass 
transition might be avoided in finite dimension be- 
cause of some still unknown mechanism. Still the 
existence of metastable glassy states seems well es- 
tablished and mean field theory provides a precise 
description of their properties. 

• Indeed, as discussed in Appendix A, in finite di- 
mension, in the limit of infinite volume and if one 
waits an infinite time, only the ideal glass state (if 
any) remains stable at finite (but arbitrarily large) 
pressure; therefore in a finite dimensional system 
stable amorphous packings exist strictly only at 
Pgcp- Here we say that a packing is stable if the 
system remains close to it for an infinite time when 
the pressure is made finite by slightly decompress- 
ing the particles. It might be useful to elucidate the 
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relation of this notion of stability with the jamming 
categories of (Doncv et at, 2007; Torquato et al, 

2000) . On the other hand, the time scale needed to 
observe the instability should diverge exponentially 
in the distance from (fccp, with an associated di- 
verging length scale (Bouchaud and Biroli, 2004); 
therefore for all practical purposes we expect the 
mean field picture to hold in finite systems. A sys- 
tematic investigation of the stability of packings as 
a function of system size would be very important 
in this respect. 

• The existence of a further phase transition has been 
recently recognized in the context of mean field 
models defined on random graphs, including for in- 
stance the hard sphere model of (Biroli and Mezard, 

2001) . It has been called freezing or rigidity tran- 
sition and is characterized by the appearance of a 
finite fraction of frozen particles in the system (Se- 
merjian, 2008; Zdeborova and Krzakala, 2007). In 
other words, at densities below the freezing tran- 
sition, even if the structure is frozen, particles can 
still diffuse out of their cages, while above the freez- 
ing transition a finite fraction of them is really stuck 
inside his cage and can only vibrate; the diffusion 
coefficient is strictly zero for these particles. This 
transition is very peculiar to mean field models: 
indeed, it has been shown (Osada, 1998) that for 
finite dimensional hard spheres the diffusion con- 
stant of a tagged particle is always strictly positive 
at finite pressure. Even if this transition is avoided 
in finite dimensional system, it might play some 
role for finite volumes and finite times. Its role in 
the context of jamming should be clarified. 

• An important result that attracted a lot of atten- 
tion in the last years is the presence of an anoma- 
lous bunch of soft modes in jammed packings at 
(pj (Majmudar et al, 2007; O'Hcrn et al, 2003; 
Silbert et al, 2005; Wyart et al, 2005a). These 
modes are related to isostaticity, they are associ- 
ated to a diverging length scale and have been re- 
lated to the Boson peak observed in glasses at low 
temperature. It would be interesting to find out 
the origin of these soft modes within the approach 
presented here. 

• Another subject of discussion in the community is 
the role of rattlers (particles that are not blocked 
by their neighbors) in jammed packings. In our 
theory the presence of rattlers is ignored, due to 
the simple form we chose for the single molecule 
density, Eq. (34). A more refined ansatz for this 
quantity could allow to compute the fraction of rat- 
tlers, and it would be interesting to compare this 
with numerical simulations. Also, in presence of 
rattlers the internal entropy of the glass should re- 
main finite at jamming (even if its derivative, the 
pressure, should diverge). It would be important 
to check this explicitely. 



• An important role might be played by 3 and more 
body interactions in the effective replicated liquid, 
at least in finite dimension. In this paper we always 
ignored these interaction. The inclusion of three 
body interaction will allow to compute the second 
order in the small cage expansion and might cure 
the unsatisfactory behavior of the Baxter resum- 
mation. We believe that this is the most important 
technical point to be studied in the future. 

• The result for the number of neighbors in d — >■ oo is 
strange and should be reconsidered with more care. 
In particular one should check whether the 1RSB 
solution is stable and look more carefully to sub- 
leading corrections, also coming from many body 
interactions. 

• Our result for the clustering transition for d —¥ oo, 
tpd ~ d/2 d , suggests that a dynamical glass transi- 
tion happens around this density. It would be very 
interesting to study the Mode-Coupling equations 
in the limit d —¥ oo to investigate this possibility. 

• The Baxter resummation can be improved in dif- 
ferent ways as discussed in the last section. In par- 
ticular there are some features of the g(r) like the 
square root singularity and the peak in r = V3D 
that we are not able to reproduce at present and 
might be captured by a more careful resummation. 

• An extension to non-spherical particles such as el- 
lipsoids should be possible and very interesting, 
since these systems show a non-trivial behavior of 
packing density with aspect ratio (Donev et al, 
2004). 

• Potentials made by an hard-core part plus a short- 
range attractive potential can be investigated: at 
the Mode-Coupling level these potentials show an 
interesting phase diagram characterized by a reen- 
trant glass transition and a glass-glass transition 
line (characterized by higher-order Mode-Coupling 
singularities). These results have been partially 
reproduced by a static HNC computation in (Vc- 
lenich et at, 2006), still with no evidence of a glass- 
glass transition. It would be interesting to see if 
better approximation schemes (such as the small 
cage approximation) could describe also the glass- 
glass transition. 

• Soft repulsive potentials such as those used in 
(Berthier and Witten, 2009a,b; O'Hcrn et al, 2002, 
2003) could be studied with this method, but this 
will require matching between the small cage ex- 
pansion for hard spheres discussed here with the 
harmonic expansion of (Mezard and Parisi, 1999a). 
This might be technically difficult. 

We hope that future work will address at least some of 
these issues. 



44 



Acknowledgments 

We wish to thank in particular the authors of (Ange- 
lani and Foffi, 2007; Berthicr and Witten, 2009b; Brito 
and Wyart, 2006; Donev et al, 2005a; van Meel et at, 
2009; Schreck and O'Hern, 2008; Skoge et al, 2006) for 
sending us their numerical data, giving the permission to 
reproduce them here and answering with patience to all 
our question about the data. 

FZ benefited a lot from the stimulating environment 
of two workshops: the "Jamming" workshop in Aspen 
Center for Physics (August 2007) and the "Dynamical 
Heterogeneities" workshop held in Lorentz Center, Lei- 
den (September 2008). He wishes to thank both centers 
for hospitality and all the participants for many useful 
discussions. 

Finally, it is a great pleasure for us to thank L. An- 
gelani, L. Bcrthier, I. Biazzo, G. Biroli, J. -P. Bouchaud, 
F. Caltagirone, A. Cavagna, P. Charbonneau, B. Coluzzi, 
O. Dauchot, A. Donev, W. Ellenbroek, G. Foffi, S. Franz, 
A. Giuliani, P. Goldbart, C. O'Hern, W. Krauth, 
F. Krzakala, J. Kurchan, J. L. Lcbowitz, A. Liu, 
R. Mari, M. Mczard, D. Ruelle, M. Schrotcr, B. Scop- 
pola, G. Semerjian, M. Tarzia, S. Torquato, P. Verroc- 
chio, P. G. Wolynes, M. Wyart, L. Zdcborova, for many 
interesting discussions that were fundamental for the de- 
velopment of this work. 



Appendix A: Metastable glassy states in finite dimension 

In this Appendix we discuss how the concepts of 
metastable glassy state and of complexity should be mod- 
ified in finite dimensional systems. The discussion will 
be brief and we refer to the original literature for more 
details. 



1. Metastable states in ferromagnetic systems 

As discussed in the body of the paper, the key feature 
of mean field theory is the existence of an exponential 
number of metastable glassy states at high density. In 
mean field these states can be defined as minima of a suit- 
able free-energy functional: the TAP functional (Mczard 
et at, 1987; Thouless et al., 1977) for spin systems, or a 
suitable density functional for particle systems (Chaud- 
huri et al, 2005; Dasgupta and Vails, 1999; Kirkpatrick 
and Wolynes, 1987a). These states are stable because 
in the infinite range case the barriers may be divergent. 
However, in finite dimensional systems, states that have 
a free energy (per particle) greater than the ground state 
can always decay towards the ground state by a nucle- 
ation process, by crossing a barrier that may diverge only 
when the energy difference with the ground state goes to 
zero. 

Therefore we should find a suitable definition of these 
states. Let us first recall the way this is done, for in- 



stance, in a ferromagnetic system. In positive magnetic 
field, the state (let us call it — ) with negative magneti- 
zation is metastable with respect to the state +. How- 
ever one can study the state — by preparing the system 
in negative field where this state is stable, and then in- 
creasing slowly the field towards its final positive value. 
After some time the state + will nucleate. Still, if one 
is far enough from the spinodal of the state — , the nu- 
cleation time is long enough and one is able to follow 
the state — into the positive field region and measure 
its properties (e.g. its magnetization). In dimension d 
the life-time of the metastable state diverges for small h 
as r(h) — exp(Ah 1 ^ d ). This example shows that once 
the correct order parameter has been identified, one can 
study metastable states by adding a suitable external 
field coupled to the order parameter in order to stabi- 
lize this state, and then change the external field follow- 
ing the evolution of the state into the region where it is 
metastable. 

Actually this can be turned into a practical computa- 
tional scheme as follows. Suppose we fix the magnetic 
field to its positive value, but then perform an expansion 
of the free energy in 1 + m, assuming that the system is 
in the metastable state — and therefore its magnetization 
to is close to —1. One can check explicitly that below the 
critical temperature for ferromagnetism the free energy 
obtained in this way displays a minimum at to = —to*, 
that gives the magnetization of the — state. One might 
wonder why a minimum appears, as the free energy for 
positive magnetic field should be convex and have a sin- 
gle minimum at positive to. The key observation is that 
the decay of the metastable state is non-perturbative in 
1 + to, so that it is missed at any order in perturbation 
theory 34 around to = — 1. Thus, the expansion around 
to = — 1 stabilizes the metastable minimum also for a fi- 
nite dimensional system. This might seem an artifact of 
the approximation, but in some sense it reflects a physi- 
cal property of the system, the existence of a metastable 
state, that is relevant for numerical simulations and ex- 
periments as far as they are able to probe this state. This 
is exactly what we are interested in. 



2. Free energy functional for metastable glassy states: 
order parameter and coupled replicas 

In the case of a 1RSB transition to an amorphous state, 
it is not possible to identify a simple order parameter 



Technically speaking the free energy is non-analytic at the point 
h = 0, being however a C°° function of h. This fact is practically 
invisible from the expansion around m = ±1. If we analytically 
continue the free-energy in the region of negative magnetic field, 
it should acquire an imaginary part proportional to r(h)~ , sig- 
naling that a sharp definition of the properties of the metastable 
phase is not possible. See Chapter 10.5 of (Parisi, 1988) for a 
more detailed discussion. 
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since the density profile of each state is amorphous and 
depends on the state. For each state, one should add to 
the Hamiltonian a specific external potential that favors 
the specific density profile of that state, but this is im- 
possible since the density profiles are not known a priori. 

A precise definition of the complexity was given 
in (Franz and Parisi, 1997; Mezard, 1999; Mezard and 
Parisi, 2000) and it is based on the same idea that we 
discussed above in the case of the ferromagnet: one cou- 
ples an external field to the order parameter in order to 
prepare the system in the metastable state. In this case 
one needs to consider a replicated system, as discussed 
in section III, since the order parameter is the overlap 
between replicas (or equivalently the cage radius). The 
derivation is the same as in section III.B.l. We introduce 
a parameter e which is coupled to the average distance 
between replicas and allows to compute the entropy as 
a function of the cage radius (see figure 7) via a Legen- 
dre transform. The physical values of A for a given e 
are the solutions of dAS(m, tp; A) = d(m — l)e and we 
are ultimately interested in the case e = 0. The function 
iS(m, ip; A) is therefore the entropy of the system of m 
copies with the constraint that each pair of copies is at 
fixed (in the thermodynamic limit) distance A given by 
Eq. (13). 

The typical shape of S(m, ip; A) is reported in Fig. 7. 
In a mean field model, at densities below ipx, the cor- 
related phase is metastable and the system is liquid; at 
the glass transition it becomes stable as discussed above. 
In mean field all this is well defined an one can perform 
explicit computations (Mezard, 1999; Monasson, 1995). 
In a finite dimensional system, however, <S(m, ip; A) must 
be a convex function of A: then the minimum at small 
A should disappear below ipx- However, if we compute 
S(m, ip; A) in a power series expansion around A = 0, 
we will find a stable minimum at small A. This is ex- 
actly the same effect that we discussed for the ferromag- 
netic case. In this way, the properties of the metastable 
state can be studied also in finite dimensional systems. 
Let us focus on the region where the equation for A 
has two solutions that correspond to a local minimum 
of S(m, ip; A) — d(m — l)eA. At fixed density, by varying 
e one of the two solutions looses its stability and it disap- 
pears: these two curves are the equivalent of the spinodal 
lines in usual first order transition. The dynamical tran- 
sition is the point where for the first time the small- A 
solution exists at e = 0: only at higher density the two 
coupled replicas may remain at an high value of the over- 
lap in absence of a force that keeps them together. On 
the contrary the static transition is characterized by the 
fact that the coexistence line of the two solutions touches 
the axis e = 0. 

In finite dimension, general arguments tell us that the 
free energy is a convex function of A, so that the correct 
shape of the function S can be obtained by the Maxwell 
construction. To discuss the consequences of this fact on 
the definition of the complexity let us consider the func- 
tion q(e) = Ao/A(e) (for some reference value of Aq) for 
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FIG. 24 The shape of the function q — Ao/A for ip > ifd- the 
full line is the correct result and the dashed line is the output 
of a mean field approximation. 

densities bigger than ipd, that is shown in figure 24. As 
can be seen from the figure, the point where we evalu- 
ate the complexity (i.e. e = and high q) is always in 
the metastable region for ip < ipx where the equilibrium 
complexity is non zero. This causes an intrinsic ambigu- 
ity in the definition of complexity because the free energy 
in not defined with infinite precision in the metastable 
phase. However we can use the fact that the free en- 
ergy is a C°° function of e near the discontinuity point 
to extrapolate the high e free energy in the metastable 
region. This ambiguity becomes smaller and smaller the 
more we approach the Kauzmann density and in general 
it is rather small unless we are very near to the cluster- 
ing phase transition. This ambiguity is not important 
from practical purposes; however it implies that there is 
no sharp, infinitely precise definition of the equilibrium 
complexity. If we forget this intrinsic ambiguity in the 
definition of the complexity we may arrive to contradic- 
tory results. 

It is important to remark once again that this discus- 
sion can be turned into a practical computational scheme 
to obtain the complexity analytically, via the small cage 
expansion, or in numerical simulations, where it has been 
used for both Lennard-Jones (Coluzzi et at, 1999) and 
hard sphere (Angelani and Foffi, 2007) systems, see the 
original papers for details. 

3. The physical meaning of the small A minimum of the 
free energy 

The previous discussion was based only on the analytic 
properties of the replicated free energy S(m,ip; A). To 
complete the description, it is desirable to have a under- 
standing of the mechanism that is beyond the metastabil- 
ity of the small A minimum. This is a key point because, 
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as we will see, the mechanism is completely different in 
mean field and in finite dimension. 



a. Mean field 

Let us first discuss what happens in mean field. In 
this case, glassy states are really stable, in the sense that 
they are local minima of a well defined free energy func- 
tional, separated by barriers whose height diverges in the 
thermodynamic limit. In a dynamical perspective, once 
the system is prepared in a glassy state, it will remain 
there forever. There is a real "ergodicity breaking" in the 
liquid state above ip^, corresponding to the ideal Mode- 
Coupling dynamical transition. So why is the small A 
minimum metastable with respect to the liquid one be- 
tween ifd and ifx? 

The reason is that the number of glassy states is expo- 
nentially large in system size between tpa and ifK- Con- 
sider for instance the case of two coupled replicas, m = 2. 
If the limit of zero coupling, if both replicas are in the 
same state, the total entropy will be twice the internal 
entropy of a state, plus the complexity X that represents 
the contribution of all possible states in which the two 
copies can be. On the other hand, if the two replicas are 
uncorrelated, the total entropy is twice the internal en- 
tropy, plus twice the complexity, since now each replica 
can be in each state independently. It is clear that by 
forcing the two replicas in the same state, one loses en- 
tropy. This is why the small A minimum has lower en- 
tropy, or higher free energy, and is metastable, as long as 
the complexity E > 0. 

Therefore, in mean field, glassy state are stable but 
their huge number is responsible for the fact that the 
replicated system finds more convenient to have each 
replica in a different state. This situation is very strange: 
indeed it is very specific of the mean field structure and it 
changes completely when one considers finite dimensional 
systems. 

b. Finite dimension 

What goes wrong in the above picture when applied to 
finite dimensional systems? 

The main problem is that, in a finite dimensional sys- 
tem, an exponential number of equilibrium states cannot 
exist. This can be argued as follows: 

1. Pure states are defined by taking the thermody- 
namic limit with a specified sequence of boundary 
conditions; and the number of different boundary 
conditions scales as e KL , i.e. as the exponential 
of the surface of the system, so that the number of 
states cannot be exponential in the volume. The 
picture of the liquid state split in an exponential 
number of stable pure states does not make sense 
in finite dimension. On the contrary, at and above 
ifiK the number of states is not exponential in the 



system size even in mean field. Therefore in this 
case the picture does not need to be modified in fi- 
nite dimensions: a sub-exponential number of pure 
states are perfectly allowed as the number of pos- 
sible boundary conditions is still very large. 

2. A more sharp argument can be done for soft par- 
ticles if we consider the minima of the energy. Of 
course there can be an exponentially large number 
of minima of the energy. A true metastable state 
would be a configuration, whose energy cannot be 
decreased by moving a finite number of particles. 
It is easy to argue that such a configuration must 
have the same free energy of the ground state. On 
the other hand we can consider local minima of 
order fc, i.e. configurations whose energy cannot 
be decreased by moving at most k particles. It is 
quite evident that the associated complexity den- 
sities will be different from zero, although 
they will tend to zero when k goes to infinity (Biroli 
and Monasson, 2000). Note that analytic computa- 
tions based on the small cage expansion discussed 
in this paper will likely predict something similar 
to Si(-E), and the effective value of k will increase 
with the precision of the computation. However for 
not too large times, and in particular in the exper- 
iments it is quite possible that Ei(-E) or T, 2 (E) are 
the relevant quantities. 

The two arguments above strongly suggest that the 
metastability of the small A minimum in finite dimen- 
sions is connected to the large scale properties of the 
states and in particular to large-scale rearrangements. 
This is at the origin of the ambiguity in the definition of 
the complexity that we discussed in the previous section. 

Then, in what sense an exponential number of 
metastable glassy states, giving a finite complexity, can 
exist in this case? A consistent "real space" formulation 
of the problem has been recently discussed in a series 
of papers (Bouchaud and Biroli, 2004; Franz, 2005; Xia 
and Wolynes, 2001a, b). One first assumes that for a fi- 
nite (and small) system it is possible to define metastable 
states and that they initially are exponentially many in 
system size. Then, one can ask whether these states re- 
main stable on increasing the system size. The key ob- 
servation is the following (Bouchaud and Biroli, 2004). 
Consider a ball of size R inside the system. The rest 
of the system will influence the ball through its bound- 
ary, therefore favoring one particular state a of the ball 
among the exponentially many, through a surface con- 
tribution F a — —kR^ 1 to the free energy of the ball 
when the latter is in state a. On the other hand, the ball 
might be in any other state (3 ^ a, and in this case it will 
gain entropy because of the large multiplicity of states: 
F not a = F(al\ possible (i ^ a) = -TY>R d . The ball will 
choose whether to stay in state a or not according to 
which is the largest between F a and F not a . It is easy 
to see that for small enough R, F a wins, while for large 
enough R, the other term wins and the ball will choose 
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not to stay in the state a suggested by its boundary. This 
reasoning shows that, for large R, the assumption of the 
existence of an exponential number of states is inconsis- 
tent, since for entropic reasons subsystems will always 
ignore the influence of their boundary and escape from 
the state. 

See (Franz, 2005; Xia and Wolynes, 2001a,b) and in 
particular (Bouchaud and Biroli, 2004) for a much more 
detailed discussion. In this series of papers the discussion 
above was precisely formulated in a nucleation theory 
where the driving force is of entropic nature. Remark- 
ably, methods to compute the length of the critical nuclei 
(droplets) and their relaxation time have been proposed; 
this is a first step towards a quantitative description of 
activated processes between ipd and ipK, even if a com- 
plete theory is still lacking. 



c. Finite dimension, finite volume 

The last point to be discussed is the assumption made 
in the previous section, that for finite and small enough 
systems, an exponential number of states still exist and 
can be taken as the starting point for the nucleation the- 
ory discussed in (Bouchaud and Biroli, 2004; Xia and 
Wolynes, 2001a). 

We will discuss this crucial point in the specific case of 
hard spheres, that is of interest here. Consider a system 
of hard spheres enclosed in a finite box with (for simplic- 
ity) hard walls. It is very easy to see that for a small box 
and high enough density, there will be disconnected sets 
of configurations, in the sense that there are pairs of con- 
figurations that cannot be transformed one into the other 
simply by moving continuously the spheres. Therefore we 
can group each set of configurations that can be contin- 
uously transformed one into each other in a "state" . We 
can say that configurations inside a state are "blocked" 
in the sense that once the system is prepared in a state it 
will never escape 35 . We can define a finite- volume com- 
plexity as £y = pTnA/y, where Nv is the number of 
such states in a volume V. 

Now we can increase the volume of the box while at 
the same time adding particles in order to keep the den- 
sity fixed. On a general ground, we expect £y to be a 
decreasing function of V since by increasing the volume 
one will open new channels to "unblock" the configura- 
tions and connect some states that will be merged in a 
bigger state. On the basis of the analysis in (Bouchaud 
and Biroli, 2004), we expect three possible behaviors of 
£y: 



We use the word "blocked" in order to avoid confusion with the 
concept of "jamming". In a jammed configuration no particle 
can move. Instead, in a blocked configuration particles can move 
a little but the whole system is unable to visit all the phase space; 
the existence of blocked configurations was shown long time ago 
by Ruelle (Ruelle, private communication). 



1. for (p < <pd, in the un-clustered liquid phase, we 
expect Xy to fall very rapidly to 0, as for large 
enough volumes all the configurations will be con- 
nected and form the unique liquid state. 

2. for ipd < f < fx we are in the clustered liquid 
phase; then we expect that £y first decreases to a 
finite plateau for small enough volume. The plateau 
correspond to the mean field value £(<p). On larger 
length scales, once the nucleation effect described 
in the previous section becomes effective, the glassy 
states are "unblocked" and the complexity will drop 
to zero; only the liquid state survives. 

3. for ip > ifK, the plateau value goes to zero since the 
mean field complexity vanishes. Therefore in this 
case the complexity will go to zero fast as in case 
1. Still, even in the V — ¥ oo limit a large number 
of amorphous states will survive (but their number 
is not exponential). 

It is interesting to remark that in case 2, "unblocking" a 
configuration will require a very large time since a large 
number of particles have to rearrange together. There- 
fore the dynamics in this region will be very slow and the 
system will stay for a long time close to a configuration 
that is made by patches of locally blocked configurations, 
the "mosaic state" of (Xia and Wolynes, 2001a). Con- 
versely, in case 3, the dynamics will be completely frozen 
even in the thermodynamic limit. 

A consequence of this discussion is that an infinite sys- 
tem of hard spheres at finite pressure, even if very large, 
will always relax to the ideal glass state, i.e. at density 
<Pgc p (if crystallization is avoided) . Therefore the ideal 
glass states at (ficcp are the only jammed states that re- 
main stable if pressure is made very large but finite, in 
the limit V — > oo. Of course for finite sizes or finite times 
this will be true also for metastable glasses. 



Appendix B: The "link" expansion for Hard Spheres 

We will describe here a formal way to justify the in- 
troduction of the effective potentials to integrate over 
m — 1 replicas and obtain an expression for the repli- 
cated entropy which is formally equal to the one of a 
non-replicated liquid. We use one replica, say replica 1, 
as a reference, and integrate over the small displacements 
of the other m— 1 replicas. Note that we are not breaking 
the replica symmetry: we are only looking for a way to 
expand the entropy for small A. In the following we will 
use the notation x = (x2, ■ ■ ■ , x m ). 

The expansion of a given diagram proceeds as follows. 
At the zeroth order (A = 0), the function p(x) is a prod- 
uct of delta functions, therefore all the replicas coincide. 
Then we have 

V = Vo = UY[pd Xl Y[(x(l)-l) , (Bl) 

J A 
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i.e. T> is the corresponding diagram of the non-replicated 
system. 

The quantities x\ — x a are of order \[A. Then the 
crucial observation is that, for a > 2, the function x( x ai~ 
x a j) in Eq. (33) is essentially constant if |ccii — aiijl differs 
from D by a quantity 3> \/A, in fact it becomes simply 
equal to xi x i — Vi) = 9{\ x i — Vi\ — D). This means 
that the integration region in the space xu such that 
\\xu — xij\ — D\ >• y/A for all links £ — (i < j) does not 
give any contribution apart from the zeroth order one. 
Let us call a link £ such that \\Xu — x\j\ — D\ <~ \J~A~ a 



The integral over 2?x contains all the vertices of the 
original diagram. However, the integrand depends only 
on the vertices that are adjacent to one of the links 
L n = (£i, ■ ■ ■ ,£ n )'i calling this set dL n , and using that 
the integrals over the other vertices give 1, we obtained 
the last equality in (B3). 

Note that, for example, the function Q2(£,£') depends 
on three or four Xi depending whether the two links £, I' 
are adjacent or not. But if the two links are non-adjacent, 
then Q 2 {£,£') = Q{£)Q{£'). This motivates the introduc- 
tion of the connected functions, defined by Q c (£) = Q(£) 
and 

Q 2 (£,£') = Q(£)Q(£') + &(£,£') , 
Q 3 (£,£',£") = Q(£)Q(£')Q(£") 

+ {Q C 2 {£,£')Q{£") + perm.} + &(£,£',£") , 

(B4) 



"critical link": the idea is to isolate the contribution of 
the regions such that k links are critical, whose volume 
is of order A k l 2 . 



1. Expansion of the replicated entropy 

We define Vx\ = \\ i pdxu and 2?x = Yli(p(i) / p)dxi 
and note that using Eq. (34) we have / Dx = 1. Defining 
Xa{£) — x( x ai — x aj), we have for a diagram V contribut- 
ing to S: 



(B2) 



(B3) 



and so on, as usual in statistical mechanics. The im- 
portant property of these functions is that they are non- 
vanishing only if the sub-diagram identified in V by the 
links in L n is connected. 

Now we insert the expression of the connected func- 
tions in (B4) in (B2). It is not difficult to check that we 
can resum the contributions containing Q{£) to obtain 

V = \ j Vx 1 ]J(xi(£) + Q(£) - l) 

O II (xi (O + -!) + ■■■ 

(B5) 

The interpretation of the above equation is the following: 
the original diagram of the replicated liquid generates a 
diagram in which the vertices carry a factor p and the 
links carry a function xi + Q — 1, plus a sum of contribu- 



V = | fX\p(i)dxud*iJl{x{£)-l) = \ j V^D Xl \{[xi{£)-l + m-Xi(£)] 

= \f v Xl J\[xi{£)-i) + \ j zteiX;iI(xi(0-i) / T>*W)-xm 

+ 4 /^iE n (xi(n -i) I v*w) xmw) - xi(n + ■■■ 

= V + fv Xl ^Q(£)H(xi(£')-l)+ fv Xl J2W,£') II (Xi(0"l) + 
•* l t'j^l Kl' t"^l,l' 

= X>0 + E / ^ E Qn(h,---,£n) I] 

„\1 J is ^ ^ I) n ifi) d \ 



having defined the functions Q n {£\, ■ ■ ■ ,£ n ) from 



Qi(£) = Q{£) - J v*[ x {£) xi W] = / dxidxj^^-lxiij) - X i(i,j)} , 
Q n (£i, ■ ■ ■ ,£ n ) = f x*f[[x(4) - xi(m = / II d ^— II - xi (£)] • 
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tions in which a subdiagram made of links L n is marked 
and on each of its connected parts an interaction Q c is 
placed. The unmarked links carry again the function 
Xi + Q ~ 1) an d the vertices a factor p. In fact, it is easy 
to check on specific examples that the explicit sum over 
the links in L n can be rearranged by grouping together, 
as usual, diagrams with the same topology. In this way 
the multiplicity factors S are corrected to become the ex- 
act symmetry factors of the diagram with marked links. 
In summary, we can write 

V = V [ Xl +Q}+V 2 [ Xl +Q, Q c 2 ]+V 3 [ Xl +Q, Q c 2 , QD+- ■ ■ , 

(B6) 

where the V n represent diagrams which are built from 
the original diagram by marking n links in all the pos- 
sible topologically inequivalent ways and placing the in- 
teractions Q c on the marked connected parts. Each di- 
agram has a symmetry factor S which is the number of 
equivalent relabeling of the vertices, taking into account 
the presence of the marked links. Clearly, summing the 
contribution of all the diagrams and Eq. (35) we get 

S[p(x),x{x,y)} = NS harm + S [p, xi{x - y) + Q(x - y)] 

+ S 2 [p, Xi{x -y) + Q(x - y), Q° 2 {x, y,z)} + --- 

(B7) 

where So[p, 1 + /] is defined in the same way as the func- 
tional (32), but for one single copy of the system and with 
/ on the links of the diagrams, while the other terms come 
from interactions involving more than two particles. 

It is worth to note at this point that the specific form 
of the effective potentials Q c n depends on the topology of 
the sub-diagram that carries them. For the three-body 
interactions there is only one possible diagrams (the one 



in which £,£' are adjacent, i.e. they share one vertex) 
but for n > 3 there are many possible diagrams that 
correspond to different interactions. 

The Q c functions are quite difficult to handle. How- 
ever, from their definition one can show that, for small 
A, they arc non-zero only if for all the links it holds 
\xi — Xj\ ~ D + 0(>/A), as expected on the basis of the 
argument put forward at the beginning of this section. 
Moreover, for m > these functions are bounded, in the 
sense that when they are different from zero, they stay 
finite for A — > 0. This implies that, for example, the 
contribution coming from diagrams containing Q c 2 is of 
order A, because the integrals over the two link variables 
have support only on an interval of order \J~A. There- 
fore, the expression (B7) is an expansion in powers of 
\[A~, where the term S n is proportional to A n l 2 . Note 
that this is specific to the hard sphere potential which is 
always constant except for the discontinuity in r = D. 
For smooth potentials, the corrections come from all val- 
ues of \xi — Xj\, and in this case it is more convenient to 
integrate over the displacement of the replicas at fixed 
center of mass, see (Mezard and Parisi, 1999a). Another 
remarkable property of the effective potentials is that in 
the limit of infinite dimension only the two body poten- 
tial should give important contributions. 



2. The effective potentials 

To gain physical intuition on the effective liquid de- 
scribed by the entropy functional (B7), it is convenient 
to compute its partition function. First of all we define 
the potentials Q by 



xi(h) ■ ■ ■ 



/II 



*f n 



IIx-W-i 



.a=2 



n 



.a=2 



(B8) 



dL„ 



and similarly for the connected potentials. The bracket indicate the average over the p(i) of the vertices in dL n . 

When constructing the partition function that generates the diagrams in (B7), we must take into account the fact 
that if a link is labeled it carries the function Q c n instead of the usual function f(£) = xi(f)[l + Q(l)\ — 1. This leads, 
using standard liquid theory, to the following grancanonical partition function: 



ZeffM = ]T z — / d N x i[ Xl (e)(i + Q(e)) [] 

N=0 ' J i KV 



1 + 



Q%{t,t') 



n 



Ki'<l" 



(i + gW)(i + Q(0) 
Q%{i,e,i") 



{l + Q{t)){l + Q{i')){l + Q{i")){ 



1 + 



(i+QW)(i+Q(f)) 



x perm 



■} 



(B9) 



where now the "links" are all possible pairs of particles. This makes it possible to identify the n body potentials of 
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the effective liquid. For instance, we have 

/ m \ 

e -4> efs {x-y) = x{x _ y){l + Q{x _ y)) = x{x _ y) I JJ x{xa _ ya) 



\a=2 



A2) 



(x-y,x-z) — \ ^ 



Ql{x - y,x - z) 



(IC=2 X(Xa - ya)x(Xg ~ Zg))^^ 



(1 + Q(x y))(l + Q(x z)) (nr= 2 X(x a y a )) XtV (UT=2 x(x a z a )) xz 



r 



(BIO) 



and so on. As before the brackets denote averages over 
the Gaussian distribution of the molecule. The entropy 
functional (B7) is obtained by expanding Z e ff in a Mayer 
series of z and then Legendre transforming to the density. 
Then one has 

S[p(x),x(x, y)] = NS harm - pNFeff , (Bll) 
where F e ff is the canonical free energy of the liquid (B9). 

3. Correlation function of the glass 

It is interesting to compute the correlation function 
of the glass state. It is the correlation function of one 
replica, see section III.B.2. It can be computed as fol- 
lows. We choose one replica and we add to the hard core 
interaction an additional potential v(x — y), only to this 
replica. Then it is easy to show that 



P 2 i n SS[p,x] 



(B12) 



v=0 



This can be done by looking at the partition function 
Z m (e), Eq. (11). If we add a potential to replica one, we 
get an additional term 

exp j-i j dxdyv(x - y)Y^S(x - x u )S(y - yij) j . 

(B13) 

Then if we differentiate with respect to v(x — y) we get 



S\nZ m (e) 



5v(x — y) 



v=0 



~pii(x,y) = ~g G {x-y) . 



(B14) 



The derivative of S is exactly the same because it is the 
Legendre transform of \nZ m (e). 

The expression of S(m, <p; A) in presence of the addi- 
tional interaction v is very simple; indeed the only modi- 
fication is xi — > Xi e ~ v -> that has to be substituted in the 
final result (B7). Note that the functions Q, Q 2 etc. do 
not contain xi explicitly because they depend only on the 
potentials of the replicas 2, • • • , m, see Eq. (B8) so they 
will not depend on v. Therefore the correlation of the 
glass is given by the correlation of the liquid described 
by the partition function (B9). 



Appendix C: The two-body effective potential 
1. Definitions 

The function Q(x, y) = Q(x — y) is defined in Eq. (B8). 
Making use of Eq. (34) we get 

Q(x -y) = 

= P~ 2 j dxdyp(x, x)p(y, y) 
= J dXdY lA (x-X) lA {y-Y) 

dtdr) lA (O lA (r))x(X+H-Y 



[J X(Xa - Va) - 1 



,a=2 



-I 



dXdY lA {x - X) lA (y -Y){q A (X- Y)^ 1 - 1} 

(CI) 



where 



<lA(r) 



d£dr]7 A (Z)~f A (r])x(r + Z 
J dr>j2A(r')x(r- r') ■ 



rj) 



(C2) 



The last equality is obtained by introducing r' = r] — £ 
and observing that as £ and rj are independent Gaus- 
sian variables with variance A, their difference is also a 
Gaussian variable of variance 2A. Following a similar 
procedure 36 , we obtain 



Q(r) = J d? l2A {r>) {q A (r- P)" 1 " 1 - l} 

We are particularly interested in the integral 
1 



(C3) 



G m (A) 



V d {D) 
1 

V d {D) 
1 

V d (D) 



dfx(r)Q(r) 
dr[q A (r) m - q A (r)} 
df[q A (r) m - X (r)] . 



(C4) 



36 We start from the last line of (CI), introduce £ 
rj = y — Y, and then r' = rj — £. 



x — X and 
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which is the contribution to the virial of the non-trivial 
part of the two body effective potential given in (BIO); 
the different expressions in (C4) are obtained starting 
either from (C3) or directly from (CI) and performing 
similar manipulations to the ones in (CI). 

We will also be interested in the following functions: 



A dG m (A) 
Fm{A)= l-m dA 
mA 1 
= l-mV d (D) 

H m{ A)^- m dG ^ A) 



j drq A (f) r 



OA 



(C5) 



dm 



m 



J drq A (r) m \nq A (r) . 



V d {D) 

In the limit nulwe have G m (A) — > while 



H m (A) -> Hi(A) 



dr In [q A {r)\ 



V d (D) 
1 

~V d (D) 



OA 

dfqA(r)\nq A (r) . 



(C6) 



In the limit m — > 0, A = am, we will show that G m {A) — > 
Go (a) which is a finite function of a. We will also show 
that 



F m {A) -+ F {a) = a 



dGo(a) 
da 



H rn (A) -> H (a) = a^^- = F (a) . 



(C7) 



Finally, we will be interested in the large d limit with the 
following scaling: 



g m (A) = lim G m (D 2 Ad- 

d—too 

T m (A) = lim F m (D 2 Ad- 

d—l-OG 



(C8) 



Unfortunately a full analytical evaluation of these in- 
tegrals is not possible. In the following we will discuss 
some particular cases in which analytical calculations are 
possible. 



(a;i,--- ,x d ). We assume, without loss of generality 
thanks to rotational invariance, that r is directed along 
and we define u = \r'\, v = \f — r'\, 
■ + x 2 the distance between r' and the 



the axis x\ in 
and R = y 1 x\ - 
axis X\. We then look to the projection r' ± of r' on the 
hyperplanc P± — (x2,-" ,x d ) perpendicular to X\\ the 
distance of r' ± from the origin in this plane is just R, 
and we can introduce polar coordinates (R, 6 d -i) on this 
plane, thus defining a set Q d -\ of d — 2 angles that spec- 
ify the position of on the sphere of radius R in P±. 
We define in this way the change to bipolar coordinates 
? -> (u,v,0d-i). 

We wish to compute the Jacobian of such transforma- 
tion. It is easy to show that 



d(u, v, 



l) 



d(X!,--- ,x d ) 



V 



Xl 


X2 


Xd 


u 


U 


u 


xi—r 


X2 


Xd 


V 


V 


V 





de d - 


1 


d(x 2 ,--- 


,Xd) 



(C9) 



The reason for the zeroes on the first column is that the 
angles Q d -\ are independent of x\ as they describe the 
position of r'j_ in P± . We can write the determinant of 
this matrix as 



0(u,M«i-i) 



d(xi,- ■■ ,x d ) 
(xi - r)R 



uv 



X2 

R 

de d 



X\R 



UV 



x d 
R 



X'l 

R 



x d 
R 



,x d ) 



d(x2,-- 

rR 

uv 



,Xd) 



X2_ 

R 



x d 
R 



d{x2,--- .x d ) 



(CIO) 



The matrix appearing in the previous equation is just the 
Jacobian matrix for the change to spherical coordinates 
in P± and its determinant is given by 



X2 x d 
R R 



d(x 2 ,--- ,x d ) 



dRd6 d - 



1 



R d - 2 J d -l(0 d -!) 



(Cll) 



From the two previous equations we obtain the desired 
Jacobian 



2. The function qA(r) and the virial coefficient G m (A) 

The function qA(r) defined in (C2) is the convolution 
of a Gaussian and the theta function x(r). This integral 
can be reduced to a one-dimensional integral by mean of 
<i-dimcnsional bipolar coordinates. This will be proved 
in next subsection: the reader who is not interested in 
technical details can safely skip this discussion. 



a. Bipolar coordinates 

To compute the convolution (C2) we introduce d- 
dimensional bipolar coordinates in the space r' — 



dp = —R^dudvJd-Ued-^dOd-! = —R^dudvdtt^ 
r r 

(C12) 

where the integral of dVL d -i is just the d — 1-dimensional 
solid angle Ct d -i. The radius R is easily expressed as a 
function of u, v by elementary geometrical considerations: 



V2u 2 v 2 + 2u 2 r 2 + 2v 2 r 2 - u 4 - v 4 - r 4 
2~r ' 



R(u, v; r) — 

' r 

(C13) 

The (positive) variables u,v,r have to respect the tri- 
angular inequalities, so the domain of integration is for 
instance u € [0, oo), v G [\r — u\,r + u]. 

Then, performing the change of variable w = v 2 — (r — 



u) 2 , and using R(u,w;r) 



w(Aru—w 



we get 
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oo rr-\-u 



UV 



q A (r) = Q d - 1 du dv— R(u, v; rf- A B{u - D) 72A (v) 

Jo J\r-u\ r 

Oroo — ± — , , ' />4ru 

d-1 



2r 



L duu ww^L 



dw 



w(4ru — w) 
4r 2 



e *z 



We now make use of the relation 



(C14) 



(C15) 



and of the definition Qd-i = ^7l^i\ to obtain the final result 



qA{r)= J D du O 



._ , ru T / ru 



(C16) 



b. Three dimensions 

Remarkably, in d = 3 the integrand simplifies and a full calculation of (^(r) is possible: 

1 



QA (r) 



erf 



du u 



D 



(r±uV 



erf 



r + D 



2 A 



+ -\ — I e 

r \ tt 



e i^ - ) +2 



r 



(C17) 



with 



erf(t) ee -H= / 
©(*) 



^[l + erf W ] = -L 



(C18) 



dxe 



Given that (JaM — > 1 for r — D ^> v^A, this allows for 
a simple numerical evaluation of the integral in (C4) to 
compute G m (A). A full computation of Q(r) is also pos- 
sible using (C3). It is also easy to check that uniformly 
in r we have 37 



q A {r) ~ 6 



D 



0{y/A) + 0(e~^) . (C19) 



c. Finite dimension: expansion in powers of \f~A~ 

In general finite dimension the integral (C16) cannot 
be explicitly evaluated. However, an expansion in pow- 
ers of \f~A similar to (C19) holds. In fact, when A is 
very small it is easy to realize that the main contribution 



37 Note that @(t) is a "smoothed" theta function. 



to qa (r) comes from the integration over the component 
of r' which is parallel to r: in the orthogonal directions 
0(\f— r'\ — D) is essentially constant and the integration 
over these components gives a correction 0(e~ D /^). 
Then we have 



Qa(t) 



f 

J — ( 



P 4A / f 

dr' . - ~ 9(r - r' — D) = Q 



D 



(C20) 

The same result can be derived from (C16) by observing 
that for A — » 0, r <~ D, one has z = ^ — > 00 and for all 
n (Abramowitz and Segun, 1965) 



e~ z VZKzI n (z) -> 1 . 



(C21) 



Thus we have, changing variables to t = ^= and s = 

u-D 
ViA 1 



QA(t) = 



ds 



D + s\JAA 
D + tVU 



1 r°° 

-= / dse- (t - s)2 = 0(t) 
V 71 " Jo 



(C22) 



The next-to-leading orders in \[~A can in principle be 
computed using the large z expansion of the Bessel func- 
tions (Abramowitz and Segun, 1965). 
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In this way we can derive the leading order contribu- 
tion to G m (A): substituting qA(t) = 9(t) in (C4) we 
get 



G m {A) = 
d-JAA 



D 



/oo 
D 



dt 



/\/iA 



D + WAA 
D 



d-1 



[6(0™ - 9(()] 



£ *[*.)--««)] (Mm), 



(C23) 



defining 



Changing variables again to t = L 7 = and s = !L 7 = and 
using 



(?) 



' 1 4- s v 4 j4 



,( a -t)V4A 



1 + 



, (C30) 



we get 



1 /"CO , — — 

~ c -tJ- dse (.-t)V^-(t-.) a 

\/7T ./n 



e t + 



.4 



(C31) 



Qo{m) 



dt[Q(t) m - G(t)] . (C24) 



One can check that the function Qo(m) is given by 
QoM-Qo(l-m) + 0((l-m) 2 ) , 

/CO 
dte(t)lne(i) , 
-oo 

close to m = 1 and that 



(C25) 



<2o(m) 




(C26) 



for m — > 0, by using the t — > — oo expansion of the error 
function. 



d. Infinite dimension 

In the limit of infinite dimension, we are interested in 
the scaling A = D 2 A/d 2 . One could then naively expect 
the small cage expansion of the previous subsection to 



work well. However, one has z 



2D 2 A 



d , and an inspec- 



tion of the large z expansion of the Bessel function I n (z) 
(Abramowitz and Segun, 1965) shows that it is indeed an 
expansion in powers of — . Then, for n = as in (C16) 
one has to resum all orders in this large z expansion. 

This can be done without difficulty, either by looking 
at the large z expansion, or by a saddle-point evaluation 
of the integral representation 

In{z) = f dtt-^e^+i) , (C27) 

and one can show that 

lim e- d2z v / 2ird 2 zl !i ^2(d 2 z) = e - & . (C28) 

Using this result in (C16), we have 



lA(r) = j D duQ) 



— — e *3 da ,„ 

e -T^r (C29) 



VIttA 



We can use this result to compute Q m (A) as defined in 
(C8). Similarly to (C23), we get 



d-i 



g m (A) = lim V4A / dt [1 + t , 

d^oo J-d/ViX V 



elt + —\ -6{t) 



AA 



AA dte 



W 4A 



-f 

J — c 



dye y 



e 



y + A 
.Vaa, 



2 

%) 



6(t) 



(C32) 



Note that if we expand the second line of (C32) for small 
A (which amounts to setting A = in the integrand) we 
obtain exactly the leading term (C23), which indicates 
that the two limits d — > oo and A — > can be safely 
exchanged. 



e. Jamming limit 

The jamming limit corresponds to m — > and A = 
am. Note that G m (A) is not analytic in A = m = 0, so 
the limits m — > and A — > cannot be exchanged. 

However, the integral (C4) is uniformly convergent in 
the jamming limit so we can exchange the limit with the 
integral. As A — >• 0, we can use the expression (C20) for 
qA(r); using then the asymptotic expansion of the error 
function for \t\ — > oo, we can show that 



lim O 



r — D\ e 



m->0,A=am y AA J ! ! 

Then wc have 



r<D, 
r>D. 

(C33) 



G (a) = lim G m (A) = — d \ drr 



rn— >0, A— am 



D 



(C34) 
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FIG. 25 The function A m (i) for different values of m com- 
pared with the asymptotic behavior for m — > and large t. 



It is also possible to show that 
dG (a) _ 



Fo(a) = a- 



da 



,. ^ dG m {A) 

0, A=am 1 — 771 



hm m — 

f0, A=am dm 



(C35) 



4a£> d 7 



Z5 



drr d ~ L (r - D) 2 e is 



In the limit d —> oo, a = D 2 a/d 2 , Eq. (C34) becomes, 
by changing variable to y — d(D — r)/D, 



Go(a) = lim G (D 2 a/d 2 ) 

D 

v 2 

dye~ y ~4z , 



lim g m (A) 

m— >0. A— am 



(C36) 



which can be also directly derived from (C32) and (C33), 
confirming that the limits d — > oo and A —> can be 
exchanged without problems. 



3. The function Q(r) in the jamming limit 

Even if qA(f) can be fully computed, the function Q[r), 
defined in (C3), cannot be computed analytically in fi- 
nite dimension. Therefore we limit ourselves here to the 
computation at first order in y/A. In this case we have 



Qa{t) ~ (^/jj Ji the integration over the components 

of r' orthogonal to f gives 1 neglecting subleading cor- 
rections and we get 



Q(r) 



dr' 



(r') 2 

e iA 



\/4ttA 
1 . ( r-D\ 



e 



r — r 1 — D 



4.4 



rn 



\\/4mA 



) 



(C37) 



where the function A m (t) is 
ds 

> 

It is easy to show that 



A m (t) = m 



e -(tVnr-«) 3 [e(«)"»- 1 _i] . (C38) 



dtA m (t) = y/m~ / ds\e(s) m - 6(s)l 

V-co (C39) 

= \fmQ {m) -> m _).o v^V4 • 



Using, for s — » — oo, 0(s) ~ 2 ^p, , , it is easy to show that 

for m -> 0, A m (0) 4l,forl«t« l/y/m, A Q (t) ~ ^j, 

and for very large t 3> 1/y/m, A m (t) ~ e~ mt . Indeed 
we have for any finite i 



/■OO 

A (i) = lim A m (t) = 2 dyye 

m-i-0 J a 



1 - V^rfe* (1 - erf(t)) 



J t»i 



1 3 

2i2 ~ 4^4 



(C40) 



The (one-dimensional) Fourier transform of Ao(t) is 
given by 



Ao(q) 



dte^A (\t\) 



l-e*f W l-erf | 



(C41) 



V^Ao 



The function A m (t) is reported for different values of m 
in Fig. 25. 



Appendix D: Derivatives of the correlation functions 

In this appendix we will show how to compute the 
derivative dg(u, v)/dlnx(%, y) that is needed to compute 
the correlation function in the first order small cage ap- 
proximation, Eq. (68). We start from the partition func- 
tion as a function of the activity 

Z[z(x),x(x,y)} = ^2 / -j^r- JJ^^II^ 1 *'^) 



N=0 



i<j 



= J exp yj dxp(x) \az(x) + - J dxdyp 2 (x,y)\nx(x,y) 

(Dl) 

where the J sign is a shorthand for J2n=oI Tvf> an< ^ 

K x ) = J2i 6 ( x - x i)i fc(x,y) = Y,^ s ( x - x i) s (y - x j)- 

From the last expression it is straightforward to show 
that 



= (K x )) = P( x ) , 



dlnZ 
<9lnz(x) 

dlaZ 1 p 2 
omx(x,y) 2 2 



(D2) 
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We define the entropy functional, the Legendre transform 
of Z: 

S[p(x),x(x,y)} = 

\nZ[z(x),x(x,y)} - J dxp{x) \nz(x) 



= max 

z(x) 



From this it follows that 
OS 



dp{x) 
dS 



In z (x) 



dlnZ _(P_ 
dlnx{x,y) dlnx{x,y) 2 



(D3) 



(D4) 



because the explicit derivative with respect to z (x) van- 
ishes due to the maximum condition. 



Therefore we need to compute 



d 2 S 



To 



simplify the notation let us write hix{x,y) = J\ 
lnx(u, v) — J 2 and lnz(x) = j(x). Then 



d 2 S d 2 \nZ 



dJ\dJ 2 dJ\dJ 2 
d 2 InZ 



dt 



d 2 lnZ dj(t) 



dJidj{t) dJ 2 



8J 1 dJ 2 



dtds 



d 2 InZ 

dJidm 



d 2 \nZ 
dj(t)dj(s) 



d 2 InZ 
dJ 2 dj(s) i 



where we computed as follows. Recalling that the 
derivatives in the last equation are done at constant p we 
can write 



d d d\nZ 

= dT 2 p{t) = 1h ~dj(t) 



d 2 InZ 

dj(t)dJ 2 



ds 



d 2 \nZ dj(s) 
dj(t)dj(s) dJ 2 



from which 



dm 

dJ 2 



ds 



d 2 \nZ 



d](t)dj(s) 



d 2 InZ 
d](s)dJ 2 



(D6) 



(D7) 



The explicit derivatives of Z that appear in Eq. (D5) 
can be explicitly related to correlation functions of the 
liquid as follows (Hansen and McDonald, 1986): 



d 2 \nZ 



dlax(x,y)dhix(u,v) 



d 2 InZ 



d 2 InZ 



dlnz(t)dlnz(s) 



d 2 In Z 



d\nz{t)d\nz{s) 



4 KP2{x,y)p 2 (u,v)) - (p 2 {x,y)) (p 2 (u,v))\ 
P 2 

— [S(x - u)8(y -v) + 8(x - v)S(y - u)]g(x, y) 

y [gi(x, y, u, v) + p- 1 [5(u - x)g 3 (y, u,v) + 3 perm.] - g(x, y)g(u, v) j 



d\nx(x 1 y)d\nz{t) 2 



[(p 2 (x,y)p(t)) ~ (h(x,y)) (p(t)>] 



y {ff 3 (x, y, t) + [p-\5{x -t) + 6(y - t)) - l}g(x, y)} 
(p(t)Ks)) - {Pit)) (p(s)) = p{5(t -s) + ph(t, s)} 



-5(t -a)- c(t,s) 



(D8) 



where the last equation follows from the Ornstein-Zernicke relation (Hansen and McDonald, 1986). 
Finally we have 

dg(u,v) 2 d 2 S l r( .. . . . . , . . 

-j] 1 \ = ~2^i 1 v^i 1 \ = o \°( x -u)o{y-v) + 6{x -v)6(y-u) }g{x, y) 

d\nx(x,y) p z omx{x,y)omx{u,v) 2 

2 

+ y {54(2, y,u, v) + p^ 1 [5{u - x)g 3 {y,u,v) + 3 perm.] - g(x,y)g(u,v)} 



4 e 

- 9 — I dtdafaxMQ + lp-iMx-Q + Siy-ty-ilgfay)} 

x {g 3 (u,v,s) + [p^{&{u - s)+5(v - s)) - l]g(u,v)} 



S(t-s) - c(t,s) 



(D9) 
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By inspection of these contribution to Eq. (68) one can 
show that only the first one is relevant to describe the 
delta peak oig(x,y), which gives Eq. (70). 



Appendix E: Scaling close to jamming 

We will show here how to compute the behavior of 
m(if,ipj) for ip — > ipj. For this we will need the asymp- 
totic behavior of Qo(m) for m — > 0; one can show that 



Tlnmp _|_ C 



(El) 



where R m and S m are C°° functions of m with R m 
I + Rim + ■ ■ ■ and S m = Sim + S^m 2 + • • • , so that 

ln<3o(™) ~ ^ln^-ilnm+y lnm + lni? m + 0(m 3/2 ) . 

(E2) 

where 0(x a ) represents a quantity which is bounded by 
Cx a for x — > 0. From this relation it follows that 



Qo(m) ~ 2m + 2 + 2 + i? m 1 ' 
d Q'o(m)_ 1 + 1 +0(m _i/ 2) 



(E3) 



dm Qo(m) 2m 2 2m 
We have then from (75): 

d m £(TO,!^) = 



2m ' 1 — m ' Qo( m ) 
d Qo(m) 



dm(l — m) 



dm Qo(m) 



= -d+0(Vm) 



On the other hand 



(E4) 



d v X{m,ip) = S'{<p)+d 
= S'(tp) + d 



d v (A*(m))-^ 2 
(A*(m))-V2 

d v [<pY(v>)] 
<pY(<p) 



(E5) 



(p Y(<p) 



{ifi - ip) . (E6) 



and we obtain, close to ipj , 

i.e. Eq. (84). 
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